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Abstract:  This  report  documents  the  development  of  the  PAR3D 
numerical  flow  model,  with  emphasis  on  modifications  incorporated  to 
facilitate  simulations  of  contaminant  dispersion  in  complex  buildings  and 
other  enclosures.  PAR3D  is  a  general-purpose  computational  fluid 
dynamics  (CFD)  code  for  predicting  three-dimensional  flow  and  transport 
in  air,  water,  and  other  incompressible  fluids.  It  includes  a  two-equation 
turbulence  model  with  adjustments  for  buoyancy,  as  well  as  transport 
equations  for  suspended  materials  (contaminants),  dissolved  gases,  and 
gas  bubbles  (in  water).  The  code  employs  parallel  processing  with 
structured  curvilinear  grids,  which  may  deform  to  accommodate  quasi¬ 
static  free-surface  displacement  in  water. 

The  CFD  applications  included  here  concern  air  flow  and  aerosol 
dispersion  inside  two  different  man-made  enclosures.  These  consist  of 
single  chambers  with  complex  internal  geometry,  but  PAR3D  is  applicable 
to  multi-chamber  enclosures  as  well.  The  reported  developments  are  part 
of  a  larger  computational  and  experimental  effort  to  characterize  the 
dispersion  of  contaminants  in  multi-room,  multi-story  buildings. 
Comparison  of  PAR3D  predictions  with  experimental  data  from  a  multi¬ 
room  test  facility  will  be  presented  in  a  separate  report. 
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1  Introduction 

Background 

The  PAR3D  code  (Bernard  et  al.  1999)  represents  a  continuation  of  work 
that  began  with  the  MAC3D  code  (Bernard  1995, 1998)  which  was  devel¬ 
oped  for  simulating  flow  and  transport  (i.e.,  dispersion)  in  deep  water. 
PAR3D  retains  all  the  capabilities  of  MAC3D,  and  it  incorporates  a  num¬ 
ber  of  improvements  and  extensions  as  well.  Foremost  among  these  is  par¬ 
allel  processing,  which  reduces  central  processing  unit  (CPU)  time  and 
increases  allowable  resolution  (problem  size)  by  using  multiple  processors 
for  executing  flow-and-transport  simulations.  Like  the  MAC3D  code, 
PAR3D  employs  multi-component  curvilinear  grids,  but  it  also  allows 
quasi-static  grid  deformation  to  account  for  free-surface  displacement  in 
water.  For  deep-water  applications  in  particular,  the  capabilities  of  PAR3D 
include  salinity  transport,  sediment  and  biomass  transport,  and  dissolved- 
gas  transfer  and  transport  induced  by  bubble  plumes  and  mechanical 
mixers. 

Water  and  air  can  both  be  regarded  as  incompressible  fluids  at  flow  veloci¬ 
ties  much  less  than  the  speed  of  sound.  Therefore,  since  PAR3D  is  formu¬ 
lated  for  incompressible  flow  in  general,  it  is  equally  applicable  for  water 
and  air  at  low  Mach  number.  With  this  fact  in  mind,  the  code  has  been 
extended  to  facilitate  simulations  of  air  flow  and  contaminant  dispersion 
in  complex  buildings  and  other  enclosures.  Internal  walls  are  represented 
as  having  either  finite  or  infinitesimal  thickness,  and  dispersed  contami¬ 
nants  can  be  reintroduced  through  supply  vents  after  being  drawn  into  the 
return  vents  of  a  heating,  ventilation,  and  air  conditioning  (HVAC)  system. 
Contaminant  release  can  be  either  continuous  or  discontinuous  (impul¬ 
sive),  and  the  boundary  conditions  for  all  variables  can  be  altered  dynami¬ 
cally  in  the  course  of  a  simulation. 

Computational  fluid  dynamics  (CFD)  is  now  a  mature  branch  of  computa¬ 
tional  physics,  which  encompasses  both  compressible  and  incompressible 
flow.  PAR3D  is  limited  to  incompressible  flow  in  particular,  which  means 
that  the  fluid  density  is  not  affected  by  pressure.  Even  so,  fluid  density 
may  be  altered  (slightly)  by  temperature,  salinity,  and  suspended  materi¬ 
als.  In  the  presence  of  gravity,  small  variations  in  density  may  create  buoy¬ 
ant  forces  that  dramatically  influence  flow  and  transport,  so  PAR3D 


ERDC/CHL  TR-07-9 


2 


employs  the  Boussinesq  approximation.  This  means  that  small  density 
variations  are  neglected  in  the  continuity  equation  but  not  in  the  momen¬ 
tum  equation.  Whatever  the  case,  the  speed  of  sound  is  infinite  for  incom¬ 
pressible  fluids,  and  this  is  the  reason  for  the  practical  restriction  to  flow 
situations  with  low  Mach  number. 

PAR3D  uses  a  finite-volume  scheme  with  structured  marker-and-cell 
(MAC)  grids  to  compute  the  transport  of  mass  and  momentum  through 
flow  regions  of  arbitrary  shape.  Structured  grids  employ  grid  cells  that  are 
numbered  with  integer  indices  (ij,k),  which  also  serve  as  the  grid-cell 
coordinates  (i,j,k)  in  a  rectangular  computational  space.  The  MAC  desig¬ 
nation  indicates  that  vector  components  are  computed  normal  to  the  cell 
faces,  whereas  scalar  quantities  are  computed  exclusively  at  the  cell  cen¬ 
ters.  Face-normal  (contravariant)  vector  components  represent  average 
values  for  the  entire  cell  face,  and  cell-centered  scalar  quantities  represent 
average  values  for  the  entire  grid  cell.  The  grid  cells  are  all  cubical  in  com¬ 
putational  space,  but  they  may  take  any  six-sided  shape  (with  planar 
faces)  in  Cartesian  (x,y,z)  space.  The  latter  feature  gives  PAR3D  much  geo¬ 
metric  flexibility  in  accommodating  flow  regions  with  non-rectangular 
boundaries. 

The  grids  are  subdivided  into  multiple  components  that  represent  rectan¬ 
gular  blocks  of  computational  {ij,k)  space,  and  each  grid  block  is  assigned 
to  a  single  processor.  Shared  information  is  communicated  between  proc¬ 
essors  via  the  standard  message-passing  interface  (MPI)  library,  which  is 
described  in  detail  by  Gropp  et  al.  (1996). 

PAR3D  uses  a  time-marching  scheme  to  solve  the  Reynolds-averaged 
Navier-Stokes  (RANS)  equations  for  incompressible  flow.  Starting  with 
potential  flow  as  an  initial  condition,  the  code  uses  discrete  time-steps  to 
advance  the  solution  toward  its  final  condition,  which  may  be  either  steady 
state  or  time-varying.  Advective  terms  in  the  RANS  momentum  equation 
are  evaluated  explicitly  (time-lagged)  to  achieve  a  better  approximation  for 
convective  transport,  and  diffusive  terms  are  evaluated  implicitly  (time- 
advanced)  to  allow  larger  time-steps.  A  Poisson  equation  for  pressure, 
obtained  by  combining  the  momentum  and  continuity  equations,  is  solved 
iteratively  during  each  time-step  to  supply  the  pressure  gradient  needed 
for  conservation  of  mass.  Transport  equations  for  other  variables  are 
solved  in  the  same  manner  as  the  momentum  equation. 
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The  influence  of  turbulence  is  incorporated  via  the  standard  K-e  turbu¬ 
lence  model  developed  by  Launder  and  Spalding  (1974).  The  model  con¬ 
sists  of  transport  equations  (which  include  production  and  dissipation 
terms)  for  the  turbulence  energy  k  and  the  turbulence-energy  dissipation 
rate  e.  These  quantities  are  used  to  calculate  eddy  viscosity  in  the  RANS 
momentum  equation,  and  eddy  diffusivity  in  the  associated  transport 
equations  for  other  variables. 

Purpose  and  scope 

This  report  documents  the  continuing  development  of  PAR3D  and  its  ten¬ 
tative  validation  for  applications  involving  particle  dispersion  in  air. 

Earlier  reports  documenting  MAC3D  emphasized  deep-water  applications 
for  reservoirs.  In  keeping  with  the  motivation  for  recent  developments, 
this  report  focuses  on  flow  and  dispersion  in  complex  (man-made)  enclo¬ 
sures.  Experimental  validation  is  limited  here  to  single-chamber  enclo¬ 
sures  with  complex  internal  geometry,  but  the  PAR3D  code  is  equally 
applicable  to  multi-chamber  enclosures.  The  current  work  was  part  of  a 
larger  computational  and  experimental  effort  to  characterize  contaminant 
dispersion  in  multi-room,  multi-story  buildings.  Comparison  of  PAR3D 
predictions  with  data  from  experiments  conducted  in  a  multi-room  HVAC 
test  facility  will  be  documented  in  a  separate  report. 

Chapter  2  outlines  the  physical  governing  equations,  and  Chapter  3  dis¬ 
cusses  their  numerical  solution.  Chapter  4  enumerates  the  various  bound¬ 
ary  conditions  that  are  imposed  on  the  flow  field,  and  Chapter  5  gives  an 
overview  of  model  operation.  Chapter  6  presents  comparisons  of  experi¬ 
mental  data  with  PAR3D  predictions  for  air  flow  and  aerosol  dispersion  in 
two  different  enclosures.  Chapter  7  offers  concluding  remarks,  and  the 
Appendix  develops  the  boundary  conditions  proposed  for  turbulence  in 
Chapter  4. 

Mode  of  technology  transfer 

As  part  of  a  larger  research  effort,  the  results  of  this  work  have  been 
applied  to  the  simulation  of  air  flow  and  contaminant  transport  in  several 
facilities,  including  an  operational  military  building.  Some  of  these  results 
will  be  disseminated  through  technical  papers  presented  at  professional 
conferences  devoted  to  the  modeling  of  fluid  flow  and  contaminant  trans¬ 
port.  Upon  publication,  this  report  will  be  made  immediately  available  to 
its  sponsor.  It  will  be  accessible  to  the  public  through  the  World  Wide  Web 
at  http://itl.erdc.usace.army.mil/library/. 
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2  Governing  Equations 

Incompressible  flow 

Subject  to  the  Boussinesq  approximation,  which  accounts  for  density 
variation  only  in  the  gravitational  acceleration  term,  the  governing 
equations  for  incompressible  flow  are  the  RANS  momentum  equation, 


Po 


d u  / — 

—  +  div(uu 
8t  v 


=  div  f  -grad  p  +  pg 


(2.1) 


and  the  continuity  equation, 


div  u  =  0 


(2.2) 


where 


p0  =  reference  density 
u  -  Reynolds-averaged  velocity  vector 
t  =  time 

f  =  Reynolds-averaged  shear-stress  tensor 
p  =  Reynolds-averaged  static  pressure 
p  =  local  density 

g  =  gravitational  acceleration  vector 


Invoking  the  widely  used  Boussinesq  hypothesis  for  turbulence,  the 
individual  shear-stress  tensor  components  zmn  take  the  form 


*mn  =M 


dum  du. 

-  +  - 


A 


dxn  dx 


m  J 


Pok#m 


(2.3) 


where 

p  =  dynamic  viscosity  =  p0  +  pr 
p0  =  molecular  viscosity 
pr  =  eddy  viscosity 

Um  and  un=  Cartesian  velocity  components  (ui,  u2,  u3)  =  (u,  v,  w) 
xm  and  xn=  Cartesian  coordinates  (xi,  x2,  x3)  =  (x,  y,  z) 
k  =  turbulence  energy  per  unit  mass 
Smn  =  unit-normal  tensor  components 
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The  symbol  Smn  represents  the  familiar  Kroniker  delta,  which  is  unity  for 
m  =  n  and  zero  otherwise.  For  practical  applications  involving  turbulent 
flow,  one  may  assume  that  jir>  >po  so  that  the  dynamic  viscosity  reduces  to 

P~Mt=Povt  (2-4) 


in  which  vr  represents  the  kinematic  eddy  viscosity.  For  transported 
quantities  (j)  other  than  momentum,  the  Reynolds-averaged  transport 
equation  is 


div 


( u  +  uip)<p]  =  div(Dipgrad<p)  +  S(<p ) 


(2.5) 


where 

u j  =  settling  velocity  (if  any)  for  <f> 


D,j,  =  eddy-diffusion  coefficient  =  — 

c j ■/,  =  Prandtl-Schmidt  number  for  (f> 
S(^)  =  source  term  for  (f> 


Turbulence 


The  k-e  model  developed  by  Launder  and  Spalding  (1974)  defines  the 
kinematic  eddy  viscosity  vr  to  be 


(2.6) 


where 

Cy  =  0.09 

k  =  turbulence  energy  per  unit  mass 
e  =  dissipation  rate  for  k 

The  coefficient  Cv  is  a  universal  (empirical)  coefficient,  and  the  governing 
equations  for  k  and  e  are  respectively  formulated  as 

—  +  div[uk)  =  P-s  +— of/V(vrgrac/k) 


(2.7) 
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5  £  £  1 

—  +  div(ue)  =  —  (c1P-c2^)  +— c//V(vT  grads 


(2.8) 


where  P  is  the  rate  of  production  for  the  turbulence  energy,  and  Ci,  c2,  ok, 
and  <je  are  empirical  coefficients. 

To  reproduce  the  observed  equilibrium  between  production,  dissipation, 
and  diffusion  for  wall-bound  shear  flow  (without  a  pressure  gradient),  as 
discussed  by  Rodi  (1980),  the  coefficients  Ci,  c2,  and  a£  must  satisfy  the 
relation, 


q  -  c2  +  - 


K 


=  0 


(2.9) 


The  symbol  k  represents  von  Karman’s  constant,  whose  value  lies  between 
0.37  and  0.44.  In  PAR3D  this  quantity  is  assumed  to  be  k  =  0.40. 

In  the  standard  k-e  model,  the  coefficients  c2,  ok,  and  oE  are  assigned  the 
values 

c2  =  1.92 

Ok  =  1.0 

Ok  =  1-3 


and  it  then  follows  from  Equation  2.9  that  the  remaining  coefficient  Ci  is 


Ci 


(2.10) 


The  rate  of  turbulence  energy  production,  which  arises  from  shear  and 
buoyancy  in  the  mean  flow,  is  given  by 


P  =  vT 


% 


2  2  2 

+  vy  +  wz 


)  +  (uy+vxf  +(uz+wxf  +(vz+wy) 


2  +  pA 


Po 


(2.11) 


in  which  the  subscripts  (x,  y,  z )  indicate  partial  derivatives  of  the  Cartesian 
velocity  components  (u,  v,  w)  and  the  fluid  density  p.  The  symbol  g 
denotes  gravitational  acceleration  (downward)  in  the  vertical  (upward) 
z-direction.  Note  here  that  stable  density  stratification  (pz<  o)  reduces 
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production,  whereas  unstable  density  stratification  (pz>  o)  has  the  opposite 
effect. 

The  k-e  model  provides  approximations  for  all  six  components  of  the 
Reynolds-averaged  shear  stress,  along  with  governing  equations  for  the 
associated  turbulence  quantities  k  and  e;  but  it  does  not  provide  boundary 
conditions  for  the  latter.  These  must  be  supplied  separately,  based  on 
other  considerations  discussed  in  the  appendix. 

Fluid  density 

The  fluid  under  consideration  (air  or  water)  is  presumed  to  be 
incompressible,  which  means  that  its  density  is  unaffected  by  pressure.  In 
spite  of  its  incompressibility,  the  effective  density  of  water  may  be  altered 
by  temperature,  salinity,  gas  bubbles,  and  suspended  materials  (sediments 
or  contaminants).  Likewise,  the  effective  density  of  air  may  be  altered  by 
temperature  and  suspended  contaminants.  For  air  and  water,  PAR3D 
imposes  a  reference  density  defined  at  a  predetermined  absolute 
temperature  T0;  i.e., 


Po  =  P{T0)  (2.12) 

For  air  at  (or  near)  atmospheric  pressure,  the  effective  density  is  inversely 
proportional  to  absolute  temperature  T  (in  degrees  Kelvin)  and  directly 
proportional  to  contaminant  mass-concentration  (f>  (in  kilograms  per  cubic 
meter)  such  that 


P~Po 


Zo 

T 


+  (p 


PqTq 
Pj  , 


(2.13) 


where 

p0  *  1.27  Kg/m3 

pc  =  bulk  density  of  contaminant  in  kilograms  per  cubic  meter 
To  =  277  degrees  Kelvin 

The  effective  density  for  water  is  likewise  proportional  to  contaminant 
concentration  (f>  and  void  ratio  9,  but  its  dependence  on  temperature  T  and 
salinity  6  is  nonlinear;  i.e., 
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p«po(l-9)  +  Pl(T)  +  p2(T,0)  +  <p  1-^  (2.14) 

V  Pc  J 

where 

p0  «  iooo  Kg/m3 
9  =  void  ratio  for  gas  bubbles 
Pi(T)  =  -0.009073 1 T-  To |  -  0.006047  O'-  To)2 
pziT, 6)  =  [1.0  -  0.001779 (T-  Ti)]  [0.80760] 

Ti  =  273  degrees  Kelvin 

The  term  “void  ratio”  refers  here  to  the  local  ratio  of  bubble  volume  to 
composite  volume  (water  and  bubbles).  The  density  variation  with  void 
ratio  and  contaminant  concentration  is  necessarily  linear,  because  gas 
bubbles  (in  water)  and  suspended  materials  (in  air  or  water)  fully  displace 
the  fluids  in  which  they  are  suspended.  In  contrast,  salt  goes  into  solution, 
and  it  does  not  fully  displace  the  water  in  which  it  is  dissolved.  Note  that 
the  salinity  0  must  be  specified  in  parts  per  thousand  (grams  per  kilogram) 
in  the  definition  for  p2(r,0).  For  temperatures  from  273  to  313  degrees 
Kelvin  (o  to  40  degrees  Celsius)  and  salinities  from  o  to  40  parts  per 
thousand,  the  combined  functions  pi(I)  and  p2(T,0)  reproduce  (within 
about  6  percent)  the  density  variation  specified  by  a  more  complicated 
equation  of  state  published  by  Fofonoff  (1985). 

Curvilinear  coordinates 

In  the  PAR3D  code,  the  governing  equations  (Equations  2.1,  2.2,  2.5,  2.7, 
and  2.8)  are  transformed  from  Cartesian  coordinates  ( x,y,z )  to  general 
curvilinear  coordinates  (£,//,£)  to  facilitate  applications  with  non- 
rectangular  boundaries.  The  one-to-one  coordinate  transformation  allows 
each  of  the  Cartesian  coordinates  to  be  expressed  directly  as  a  function  of 
the  curvilinear  coordinates,  such  that 
x  =  x  (f,i7,£) 
y  =  y  (f,i i,0 
z  =  z  (£17,$ 

It  then  follows  that  the  x-,  y-,  and  z-derivatives  of  any  variable  <f)  are 
related  to  its  y-,  and  ^-derivatives  by  the  chain  rule.  Thus,  the 
transformed  x-,  y-,  and  z-components  of  grad  ^become 


(2.15) 
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<Py=Zy<Pz+Vy<Pv+Cy<P( 

(2.16) 

<Pz=tz<P£+1h(P,1+Cz(Pc 

(2.17) 

Here  subscripts  ( x,y,z )  and  (£47,  £)  indicate  derivatives  with  respect  to  the 
Cartesian  and  curvilinear  coordinates,  respectively.  As  discussed  by 
Anderson  et  al.  (1984),  the  Cartesian  derivatives  of  the  curvilinear 
coordinates  are  given  by 


4=^(y^-y^) 

(2.18) 

4  =  -J1vr¥,) 

(2.19) 

(2.20) 

= -J 1  (y  tzc  -  y  cz^ ) 

(2.21) 

%=J-1(xfr~xczs) 

(2.22) 

Tiz  =  -J~1(xfye-xeys) 

(2.23) 

(2.24) 

Cy=-J1{xizv-xrizi) 

(2.25) 

C  =  J_1(x#y,-v^) 

(2.26) 

The  symbol  J  represents  the  Jacobian  determinant  of  the  coordinate 
transformation, 
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Expressed  in  expanded  form,  J  becomes 

J = xf  (y,z<  -  yczn )  -  \  (y,z;  -  y& ) + xc  ( y *z,  -  Vf )  (2-28) 

Under  the  coordinate  transformation  outlined  above,  the  curvilinear 
expression  for  the  divergence  of  any  Cartesian  vector  u  =  (u,v,w )  is  given 
by 


div  0  =  U1  ( 

where 

u  =  j(4u- 
V  =  J(?]XU- 
w  =  j(Cxu- 


f.+^+W,) 

(2.29) 

4yV  +  ZzW) 

(2.30) 

rjyv  +  r]zw) 

(2.31) 

Cyv  +  Czw) 

(2.32) 

In  PAR3D  the  Cartesian  velocity  components  ( u,v,w )  are  retained  (for 
convenience)  as  dependent  variables  in  the  transformed  RANS 
momentum  equation,  but  the  transformed  velocity  components  (w,v,  w) 

are  used  in  the  transformed  continuity  equation.  The  latter  represent  the 
proper  components  for  the  time-dependent  velocity  field  in  the  curvilinear 
coordinate  system,  but  their  time  increments  are  constructed  from  the 
more  conveniently  computed  increments  of  the  Cartesian  components,  as 
discussed  in  Chapter  3. 
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3  Numerical  Considerations 

Explicit  upwind  scheme  for  advection 

All  transport  equations  in  the  PAR3D  code  retain  the  conservative  form 
div(u  <f)  for  advective  (i.e.,  convective)  terms.  In  principle  the  continuity 
equation  (Equation  2.2)  reduces  these  to  u  •  grad  <fr,  but  in  practice  the 
conservative  form  achieves  a  better  representation  of  advective  transport 
on  non-uniform  grids.  In  any  case,  a  one-dimensional  (l-D)  example  is 
used  here  for  illustrating  the  conservative  discretizaton  of  advective  terms 
in  general.  Consider  the  following  l-D  equation  for  transport  by  advection: 

=  0  (3.1) 

at  orj 


where 

(f>  =  transported  quantity 
t  =  time 

v  =  transporting  velocity  in  the  77-direction 

Let  the  time  and  space  increments  be  At  and  A77,  respectively;  and  let  the 
integer  subscript./  indicate  discrete  values  (f>j  defined  at  stations  j-i,j,j+i, 
etc.,  along  the  77-axis.  Furthermore,  let  the  integer  superscript  n  indicate 
discrete  values  of  ^at  time  levels  n-i,  n,  n+i,  etc.  Imposing  a  forward 
approximation  for  the  time  derivative  of  (/>,  along  with  a  central 
approximation  for  the  space  derivative  of  V(j),  one  obtains  a  discrete 
approximation  for  Equation  3.1  at  time  level  n  and  position  j, 

|  °j+l/2(Pj+l/2  ~  °j-l/2(Pj-l/2  _  q  ^ 

At  A77 

If  position  j  is  regarded  as  the  center  of  a  l-D  grid  cell,  with  positions 
j  ±  1/2  as  the  corresponding  cell  faces,  then  Equation  3.2  represents  a 
finite-volume  approximation  for  Equation  3.1  in  l-D  space. 
Inconveniently,  however,  the  dependent  variable  (f>  has  been  discretized 
only  at  integer  positions  j,  and  its  intermediate  (cell-face)  values  at 
positions  j  ±  1/2  can  be  obtained  only  by  interpolation  or  extrapolation. 
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PAR3D  obtains  these  values  by  second-order  upwind  extrapolation.  For 
example,  when  t>"+1/2  >  0  then 


<1/2 =!?;-§<!  o-s) 

But  when  l>"+1/2  <  0  then 

<^§<1  (3-4) 

In  the  parlance  of  numerical  analysis,  this  is  an  explicit  upwind  scheme, 
which  is  first-order  correct  in  time,  second-order  correct  in  space,  and 
marginally  stable  for  C  <  2/3,  where  C  is  the  Courant  number, 


PAR3D  uses  second-order  upwind  extrapolation  (Equations  3.3  and  3.4) 
for  all  dependent  variables  except  k  and  e,  for  which  it  uses  first-order 
upwind  extrapolation  instead.  The  nonlinear  production  term  P  renders 
Equations  2.7  and  2.8  more  prone  to  instability  than  the  other  governing 
equations,  but  first-order  extrapolation  keeps  their  discrete  solutions 
numerically  stable. 

Implicit  central  scheme  for  diffusion 

When  expressed  in  expanded  form,  the  shear-stress  terms  arising  from 
div  f  in  Equation  2.1,  and  likewise  the  diffusion  terms  in  Equations  2.5, 
2.7,  and  2.8,  all  include  terms  of  the  form  D $  div  grad  (/).  These  terms  are 
dissipative,  and  they  contribute  to  numerical  stability  when  C  <<  1;  but 
they  can  also  hamper  computational  speed  by  reducing  the  allowable  size 
of  the  time-step  At.  Consider,  for  example,  the  l-D  advection-diffusion 
equation  (ADE): 


d(p  dm  _  d2m 

—  +  v—  =  D — V 
dt  dg  dg 


(3.6) 


Here  the  coefficient  D  represents  the  diffusion  coefficient,  and  the  other 
quantities  are  the  same  as  in  Equation  3.1.  Using  a  three-point  upwind 
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approximation  (assuming  v  >  o)  for  the  advective  term,  and  a  central 
approximation  for  the  diffusion  term,  the  explicit  approximation  for 
Equation  3.6  becomes 


Atf  +v 

'3<Pj-^<Pj.x  +  <PU~ 

=  D 

"<,-2  <Pj+<pU 

At  2 

Arj 

(A??)2 

As  a  rough  rule  of  thumb,  one  can  expect  Equation  3.7  to  be  numerically 
stable  as  long  as 


-C  +  2r<l 
2 


where  C  is  given  by  Equation  3.5,  and 


DM 

(A/7)2 


(3.8) 


(3.9) 


The  ADE  stability  rule  (Equation  3.8)  depends  on  both  C  and  T,  but  the 
advective  component  is  inversely  proportional  to  At],  whereas  the  diffusive 
component  is  inversely  proportional  to  (A77)2.  This  means  that  the  time- 
step  allowed  by  the  advective  term  decreases  with  the  grid  spacing,  but  the 
time-step  allowed  by  the  diffusive  term  decreases  with  the  square  of  the 
grid  spacing.  For  turbulent  flow  simulations,  where  diffusion  can 
dominate  advection  (locally)  by  an  order  of  magnitude  or  more,  the 
diffusive  time-step  limit  may  be  far  more  confining  than  the  advective 
(Courant)  limit.  While  it  may  be  desirable  to  use  an  explicit  (time-lagged) 
approximation  for  advection  terms,  to  achieve  better  accuracy  for 
advective  transport,  it  is  likewise  expedient  to  include  an  implicit  (time- 
advanced)  contribution  from  diffusion  terms,  to  circumvent  the  diffusive 
time-step  limitation.  The  discretization  employed  by  PAR3D  is  equivalent 
to  the  following  two-level  representation  for  the  ADE: 


^  0 

1 

00 

-<3 
-  ^ 

1 

73 

7  3 

+ 

-<3 

7  = 

ro 

_ 1 

=  D 

"<1 -2^;+1+Ci " 

At  2 

At] 

(A/7)2 

where  <7" 11  =  <7"  +  A <p" . 
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Since  implicit  treatment  of  the  diffusion  term  introduces  A (pf  on  both 
sides  of  Equation  3.10,  it  follows  that  multiple  iterations  are  required  for 
A <f>jn  in  each  time-step.  For  this  purpose,  PAR3D  uses  a  symmetric  Gauss- 
Seidel  scheme,  each  iteration  of  which  updates  A </>jn  (cell-by-cell)  in  a 
forward  sweep  from  the  origin  to  the  uppermost  corner  of  the  grid, 
followed  by  a  reverse  sweep  from  the  uppermost  corner  back  to  the  origin. 
The  aim  of  this  process  is  numerical  stability  rather  than  iterative 
convergence  for  A <f>jn  at  time  level  n.  Therefore,  only  a  few  Gauss-Seidel 
iterations  are  needed  within  each  time-step. 

Velocity  interpolation  on  MAC  grids 

PAR3D  uses  MAC  grids  that  are  staggered  in  their  discrete  placement  of 
vector  and  scalar  quantities.  This  means  that  the  discrete  values  for 
variables  such  as  pressure  and  temperature  are  formally  defined  only  at 
grid-cell  centers,  whereas  the  normal  (contravariant)  components  of 
vectors  such  as  velocity  are  formally  defined  only  on  the  cell  faces  from 
which  they  emanate.  The  most  convenient  locations  for  computing  velocity 
increments  are  the  cell  centers;  but  this  requires  an  additional  closure 
relation  that  defines  cell-centered  velocity  components  in  terms  of  face- 
centered  velocity  components,  and  face-centered  velocity  increments  in 
terms  of  cell-centered  velocity  increments.  The  following  l-D  example 
illustrates  the  closure  relation  that  is  used  in  PAR3D. 

Let  the  integer  locations  j-i,j,j+i  represent  cell  centers  along  the  ^-axis, 
and  let  the  half-integer  locations  j  ±  1/2  represent  cell  faces.  The  discrete 
velocities  vnj±  1/2  are  formally  defined  only  on  cell  faces;  but  when  they  are 
needed  at  the  cell  centers,  they  are  computed  by  simple  averaging;  i.e., 

u;"a|t<v2+<v=]  I3-11* 

The  velocity  increments  Avf  are  formally  defined  only  at  the  cell  centers; 
but  when  they  are  needed  to  update  face-centered  velocities,  they  are 
projected  onto  the  cell  faces  by  fourth-order  interpolation;  e.g., 

A  <£1/2  =  ^  [9AO?,,  -  AU;,2  +  9  A  u;  -  AU;_,  ]  0.12) 
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The  interpolated  increments  defined  by  Equation  3.12  are  the  ones  that  are 
used  to  advance  the  face-centered  velocities  from  one  time  level  to  the 
next;  i.e, 


C/2  =  C/2  +£[9AC -aC  +9A«; -4^]  (3.13) 

The  closure  relation  (Equation  3.12)  between  face-centered  and  cell- 
centered  velocity  increments  retains  the  convenience  of  cell-centered  time- 
discretization  for  all  variables.  That  is,  the  time-dependent  governing 
equations  (Equations  2.1,  2.5,  2.7,  and  2.8)  are  solved  numerically  as 
though  they  were  all  discretized  only  at  the  cell  centers.  In  the  PAR3D 
analog  for  Equation  3.12,  increment  Avnj+i/2  gives  way  to  the  face-centered 
curvilinear  increments  ( Au,Av,Aib )  computed  from  the  cell-centered 

Cartesian  increments  (Au,  An,  Ain).  In  the  transformed  continuity 
equation,  curvilinear  velocity  components  (u,n,ih)  are  formally  defined 

only  on  the  cell  faces,  and  these  are  the  velocity  components  that  represent 
the  evolving  flow  field.  In  any  case,  the  formal  definition  of  velocity 
components  only  on  cell  faces  precludes  numerical  oscillations  that  may 
arise  from  the  continuity  equation  when  velocity  components  are  formally 
defined  at  the  cell  centers  instead. 

Application  of  the  pressure  gradient 

In  the  manner  discussed  in  previous  reports  (Bernard  1995, 1998)  the 
discrete  pressure  gradient  is  omitted  from  the  discretized  momentum 
equation,  only  to  be  computed  and  imposed  after  provisional  (cell- 
centered)  velocity  increments  have  been  projected  onto  the  cell  faces.  This 
pressure  gradient  enforces  mass  conservation,  and  it  lags  the  velocity  by 
one  time-step.  For  illustration,  consider  the  inviscid  Euler  equation  with 
unit  density, 


—  +  div{uu)  =  -grad  p  (3.14) 

dt 

and  let  u  be  constrained  by  the  continuity  equation  (Equation  2.2).  The 
existing  velocity  un  at  time  level  n  is  used  first  to  compute  a  provisionally 
updated  velocity  u* ,  given  by 


u*  =un  -At div{unun) 


(3.15) 
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where  div  u*  ^  0 

The  updated  (mass-conserving)  velocity  un+1  is  then  computed  as 

un+1  =  u*  -  At  grad  pn  (3.16) 


where  the  constraint  div  un+1  =  0  requires  thatp"  satisfies  a  discrete 
Poisson  equation  for  pressure, 


„  div  u 

div  grad  p  = -  (3.17) 

At 

PAR3D  employs  a  conjugate-gradient  iteration  scheme  developed  by 
Kapitza  and  Eppel  (1987)  to  solve  the  Poisson  equation  forpn.  The  same 
iteration  scheme  is  used  to  solve  a  similar  equation, 

div  grad  ®  =  div  u*  (3.18) 

where  0*  now  represents  some  arbitrarily  specified  initial  velocity,  and  the 
required  mass-conserving  initial  velocity  for  time  level  n  =  o  is  given  by 

u°  =  u*  -  grad  <t>  (3.19) 


Computational  geometry 

To  facilitate  the  implementation  of  curvilinear  coordinates  with  structured 
MAC  grids,  PAR3D  identifies  the  coordinates  (£,//,£)  with  the  integer 
indices  ( i,j,k )  such  that 


£  =>  i 
rj  =>  j 
C  =>  k 

These  grids  may  have  almost  any  shape  in  Cartesian  space,  with  arbitrary 
spacing  (Ax,Ay,Az)  between  the  grid  nodes.  Regardless  of  shape,  however, 
they  have  unit  spacing  in  computational  space;  i.e., 


Ag  =  Ai  =  1 
A  77  =  Aj  =  1 
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=  Ak  =  1 


Figure  l  shows  a  single  non-rectangular  grid  cell  in  Cartesian  space  and 
the  same  cell  in  computational  space,  where  it  is  perfectly  cubical.  To 
further  illustrate  the  correspondence  between  Cartesian  and 
computational  spaces,  Figure  2  shows  a  grid  surface  with  non-rectangular 
edges  in  the  Cartesian  (x,  y)  plane  and  the  same  grid  surface  projected 
onto  the  computational  (i,j)  plane. 


Cartesian  (x,  y,  z)  coordinates  Computational  (/,  j,  k )  coordinates 


Figure  1.  Grid  cells  in  Cartesian  and  computational  spaces. 

Grids  may  be  subdivided  into  two  or  more  contiguous  components,  which 
may  have  any  shape  in  Cartesian  space,  but  which  have  rectangular  shapes 
in  computational  space.  These  components  are  called  grid  blocks,  and 
grids  containing  more  than  one  block  are  called  multi -block  grids.  Each 
block  occupies  a  distinct  region  of  Cartesian  space,  which  resides  on  a 
separate  computer  processor  with  its  own  (i,j,  k)  coordinate  system. 
Figure  3  shows  two  (contiguous)  grid  blocks  with  cutouts  that  represent 
internal  solid  obstacles. 


8 
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BLOCK  1  BLOCK  2 
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Figure  3.  Two  contiguous  grid  blocks  with  cutouts  representing  solid  obstacles. 


Finite-volume  discretization 

As  a  result  of  the  unit  grid  spacing  in  computational  space,  the  volume  of 
each  grid  cell  in  Cartesian  space  is  equal  to  the  value  of  the  Jacobian  J  for 
the  coordinate  transformation,  which  is  given  by  Equation  2.28.  For  each 
cell  face,  one  can  construct  an  area  vector  Ssm  normal  to  that  face.  This 

vector  is  the  product  of  the  (scalar)  cell-face  area  with  a  unit  vector  normal 
to  the  face,  and  the  subscripts  m  =  1,  2,  3  correspond  to  the  77-,  and 
£-faces,  respectively.  The  Cartesian  components  of  the  area  vectors  are 
Ssmn,  in  which  the  subscripts  n  =  1,  2,  3  correspond  to  the  x-,  y-,  and 
z-directions,  respectively.  Using  condensed  notation,  the  definition  for 
these  components  is 


Ss 


mn 


(3.20) 


where  (&,  f2,  £3)  =  (£  V>  0  and  (*i.  *2,  x3)  =  (x,  y,  z). 


For  example,  the  x-component  for  area  vector  (normal  to  a  face  of 
constant  £)  is 


te11=J  —  =  ynz(-ycz,j  (3.21) 

Scalar  quantities  (pressure,  temperature,  density,  etc.)  are  formally 
evaluated  at  the  cell  centers,  but  their  values  actually  represent  cell- 
averaged  values.  Vector  fluxes  (mass  flux,  momentum  flux,  etc.)  represent 
face-integrated  quantities  evaluated  on  each  cell  face.  Subject  to  the  finite- 
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volume  discretization  outlined  above,  the  components  of  volumetric  flux 
normal  to  faces  of  constant  q,  £  respectively,  are 


u  =  Ss±1  u  +  Ss12  v  +  Ss13  w 

(3.22) 

v  =  Ss21  u  +  Ss22  v  +  Ss23  w 

(3.23) 

w  =  Ss31  u  +  Ss32  v  +  Ss33  w 

(3.24) 

Note  that  Equations  3.22,  3.23,  and  3.24  are  equivalent  to  Equations  2.30, 
2.31,  and  2.32,  respectively.  With  u  =  0 ,  the  curvilinear  (finite-volume) 

form  of  the  transport  equation  (Equation  2.5)  now  becomes 

J~ft+Jdiv(u<p)  =  J  div(D<p  &rad  <p)  +  -1  S (<p) 

(3.25) 

The  curvilinear  form  for  the  advective  term  is 

Jdiv(u<p)  =  ^(u<p)  +  -^(v<p)  +  -^(w<p) 

(3.26) 

For  a  grid  cell  centered  at  (£,  q,  £)  =  (i,j,  k),  with  unit  volume  in 
computational  space,  the  evaluation  of  Equation  3.26  is  straightforward; 
i.e., 

A(up)  =  (u^)M/2-(u^)M/2 

(3.27) 

J,m=mW2-mW2 

(3.28) 

8 

— ( W(p)  =  (wrp )  -(w<p) 

\  '  /  k+i/2  V  '  'k-1/2 

(3.29) 

The  diffusion  term  in  Equation  3.25  is  evaluated  in  expanded  form;  i.e., 

J  div^Dp  grad  p)  =  J  grad  Dp  •  grad  p  +  J  Dp  div  grad  (p  (3.30) 
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The  calculation  of  the  dot  product  in  Equation  3.30  is  straightforward, 
using  Equations  2.15  to  2.17.  The  second  term  is  evaluated  as  follows: 


JDip  div(grad(p)  =  Dv> 


'dPi  [  dq>2 

v  dg  dr/ 


arj 


(3.31) 


where 


Ci  (Px  +  dsi2 (Py  +  c>s13 <7Z  (3.32) 

dp2  =Ss21(px  +Ss22(py  +Ss23(pz  (3.33) 

^3  —  ^31  +  <5s32^y  +  <5s33^>z  (3.34) 

and  Equations  2.15  to  2.17  are  used  to  compute  the  x-,  y-,  and  z-derivatives 
in  Equations  3.29  to  3.31.  Bernard  and  Kapitza  (1992)  explain  in  detail  the 
finite-difference  approximations  needed  for  the  77-,  and  ^-derivatives 
that  occur  in  grad  ^  via  Equations  2.15  to  2.17  and  Equations  3.32  to  3.34. 
Suffice  it  here  to  say  that,  in  these  equations,  the  computed  derivatives 
must  employ  certain  combinations  of  (interpolated)  face-centered  and 
cell-centered  coordinates  (x,  y,  z)  and  dependent  variables  ^to  ensure  that 

curl  grad  (p  =  0  (3.35) 

The  foregoing  comments  and  constraints  concerning  the  div  and  grad 
operators  apply  wherever  these  operators  occur  in  the  governing  equations 
outlined  in  Chapter  2. 
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4  Boundary  Conditions 

Solid  boundaries 

In  the  context  of  PAR3D  flow  simulation,  any  solid  boundary  that  fully 
confines  the  flow  is  considered  impermeable.  It  might  be  a  wall,  floor,  or 
ceiling  of  a  room  in  a  building;  or  it  might  be  a  bottom  or  sidewall  of  a 
reservoir  or  an  open  channel.  Nothing  except  heat  can  pass  through  a 
boundary  of  this  kind,  and  the  component  of  velocity  normal  to  the 
boundary  is  identically  zero.  In  computational  (i,j,  k)  space,  where 
(£,  r],  £)  =  {i,j,  k),  the  normal -velocity  constraint  is  easily  stated  as 

u  s  o  for  i-boundaries  (constant  £,  constant  i ) 

v  =  o  for  ^-boundaries  (constant  77,  constant/) 

w  =  o  for  k-boundaries  (constant  £  constant  k ) 

Subject  to  the  constraint  that  nothing  can  pass  through  the  boundary,  the 
gradient  flux  of  any  transported  quantity  (f>  (other  than  heat)  must  also  be 
zero  normal  to  the  boundary,  such  that 

r\*gradcp  =  0  (4.1) 

where  n  is  a  vector  normal  to  the  boundary,  as  defined  by  Equation  3.20. 
Note  that  the  foregoing  constraints  apply  to  slip  (frictionless)  and  no-slip 
(frictional)  boundaries  alike.  In  the  case  of  heat  flux,  however,  which  is 
assumed  to  be  proportional  to  the  gradient  of  temperature  T  in  the  (solid) 
material  adjacent  to  the  boundary,  Equation  4.1  gives  way  to 


q  =  -CT  n»gradT  (4.2) 

where  q  is  the  heat  flux  through  the  boundary  ( into  the  fluid)  and  Cris  the 
thermal  conductivity  for  the  material  adjacent  to  the  boundary. 


In  principle,  a  slip  wall  should  impose  no  shear  stress  on  the  flow  adjacent 
to  the  wall.  In  PAR3D  this  condition  is  idealized  with  the  tangential- 
velocity  constraint, 
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n»gradu  =  0 


(4.3) 


where  v  represents  here  the  component  of  velocity  tangent  to  the  wall.  For 
a  no-slip  wall  adjacent  to  laminar  flow  in  particular,  PAR3D  approximates 
the  shear  stress  tw  tangent  to  the  wall  as 


r 


w 


Mo n 

Ini 


•  grad  v 


(4.4) 


Subject  to  the  no-slip  condition,  which  imposes  zero  velocity  at  the  wall, 
Equation  4.4  gives  way  to  the  first-order  discrete  approximation, 


T 


w 


2  Mo  ^1 

~x~ 


(4.5) 


where 

/Li0  =  molecular  viscosity 

Vi  =  average  tangential  velocity  in  grid  cell  adjacent  to  wall 
<5i  =  thickness  of  adjacent  grid  cell  normal  to  wall 

For  a  no-slip  wall  adjacent  to  turbulent  flow,  the  shear  stress  at  the  wall  is 
proportional  to  the  square  of  the  friction  velocity;  i.e., 

rw  =  p0  vt  (4.6) 

For  “hydraulically  rough”  walls,  the  friction  velocity  v*  can  be  computed 
directly  from  the  roughness  height  6*  of  the  wall  surface  material.  PAR3D 
uses  the  following  profile  for  v  near  the  wall,  based  on  experimental  data 
cited  by  Kuethe  and  Schetzer  (1967): 

—  ~8+  for  S+<  5.25  (4.7a) 

n* 


-*-  +  0.524^  for  5.25  <S+  <20  (4.7b) 

n*  k 

—  «  — ln(9<5v)  for  £+>20 

O*  K 


(4.7c) 
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in  which  von  Karman’s  constant  is  taken  to  be  k  =  0.4. 


For  rough  walls,  the  non-dimensional  distance  8+  is  obtained  from  the 
roughness  height  8*  and  the  normal  distance  8  from  the  wall,  such  that 


10  S 

3 ~K 


(4.8) 


For  “hydraulically  smooth”  walls,  however,  the  friction  velocity  v*  must  be 
computed  from  Equation  4.7  (a,  b,  and  c)  using  molecular  viscosity  y0  with 
an  alternate  definition  for  8+;  i.e., 


e,  Po  5 

°+=  - 

Po 


(4.9) 


Also  needed  for  turbulent  flow  are  boundary  conditions  for  the  turbulence 
energy  k  and  its  dissipation  rate  e.  As  discussed  by  Rodi  (1980),  near-wall 
values  for  these  variables  (adjacent  to  no-slip  boundaries)  are  commonly 
obtained  from  the  “wall”  functions, 

for  *^1  (4.10) 

k v 

3 

s  =  —  for  S  ~S1  (4.11) 

K  l\ 

Equations  4.10  and  4.11  do  not  represent  true  boundary  conditions 
because  they  are  not  imposed  at  the  wall  itself.  They  might  better  be 
described  as  off-the-wall  surrogates  for  near-wall  modeling  with  authentic 
(properly  formulated)  boundary  conditions.  In  principle,  they  are 
applicable  only  for  fully  developed,  unidirectional  shear  flow  without  a 
pressure  gradient;  but  they  are  typically  employed  without  regard  for  flow 
geometry  or  pressure  variation.  When  used  with  coarse  grids,  which 
cannot  resolve  the  near-wall  velocity  profile,  Equations  4.10  and  4.11  may 
even  contaminate  the  numerical  solution  with  inappropriate  values  for  k 
and  e  in  the  flow  field  away  from  the  wall. 

Instead  of  wall  functions  imposed  in  the  flow  adjacent  to  no-slip 
boundaries,  PAR3D  uses  boundary  conditions  that  reproduce  the  off-the- 
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wall  values  given  by  Equations  4.10  and  4.11,  but  only  for  the 
circumstances  under  which  they  are  valid.  Specifically,  in  the  wall-adjacent 
grid  cells,  the  near-wall  turbulence  production  rate  P  is  approximated  as 

P«P„+|4-  (4.12) 

2  ko1 

in  which  P0  is  the  production  rate  otherwise  given  by  Equation  2.11.  The 
flux  components  for  k  and  e  normal  to  the  wall  are  respectively 
approximated  as 

- •  grad  k  «  — ^ — |n|(k-k*)  (4.13) 

°k  Cv  (Tk 

vTn  k  u  k  ,  , ,  . 

-*— • grads  — n  (f-f*)  (4.14) 

<j  C  <j 

s  v  s 

where 

k  =  0.40 
Cv  =  0.09 

Ok  =  1.0 
Oe  =  1-3 


7  a3 

£*  =  - 

kS j 


The  near-wall  production  rate  and  the  wall-normal  flux  components  for  k 
and  e  were  developed  in  an  unpublished  supplementary  effort,  which 
involved  the  l-D  numerical  solution  of  the  RANS  momentum  equation 
(Equation  2.1)  and  the  turbulence  model  equations  (Equations  2.7  and 
2.8)  for  unidirectional  shear  flow  without  a  pressure  gradient,  parallel  to  a 
no-slip  wall.  These  developments  are  discussed  in  detail  in  the  appendix. 

When  supplemented  with  Equations  4.6,  4.7,  and  4.12  to  4.14,  the 
numerical  solutions  for  Equations  2.1,  2.7,  and  2.8  exhibit  decreasing  wall 
influence  with  increasing  wall-normal  grid  spacing  §1.  Moreover,  as  v* 
approaches  zero,  Equations  4.13  and  4.14  reduce  to  null  conditions  for  the 
wall-normal  flux  of  k  and  e.  This  effect  is  very  different  from  that  of 
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Equations  4.10  and  4.11,  which  impose  null  values  for  k  and  e  as  v* 
approaches  zero. 

PAR3D  designates  all  grid-block  interfaces  and  internal  cutout  surfaces 
(Figure  3)  as  slip  boundaries  by  default.  Model  users  may  designate  these 
boundaries  individually  as  slip  or  no-slip,  and  they  may  also  change  the 
default  to  no-slip,  via  the  PAR3D  input  files.  To  facilitate  coarse-grid 
simulations  for  complex  enclosures  with  many  internal  partitions,  the 
code  also  allows  any  grid  surface  to  be  designated  as  an  impermeable  slip 
boundary.  In  addition,  the  grid-block  interfaces  may  be  individually 
designated  as  fully  permeable  connecting  surfaces  (called  “grid  cuts”) 
between  contiguous  flow  regions. 

Inflow  and  outflow  boundaries 

Boundaries  through  which  flow  may  enter  or  leave  the  computational  grid 
are  designated  either  as  flux  boundaries  or  as  open  boundaries.  These  may 
consist  of  groups  of  one  or  more  cell  faces  on  any  one  of  the  six  bounding 
surfaces  for  a  particular  grid  block  or  internal  cutout  (Figure  3).  On  a  flux 
boundary,  the  total  rate  of  flow  and  the  distribution  of  flow  are  fixed  for 
the  duration  of  the  simulated  time-interval.  On  an  open  boundary,  the 
total  rate  of  flow  is  fixed,  but  the  distribution  of  flow  varies  with  conditions 
just  upstream  of  the  boundary  in  question. 

Phantom  grid  cells  that  lie  outside  the  grid  or  inside  internal  cutouts  are 
called  “out”  cells.  Every  grid  block  is  completely  surrounded  by  out  cells 
that  may  contain  either  updated  values  for  adjacent  variables  in 
neighboring  grid  blocks  or  updated  values  for  variables  entering  the  grid 
via  local  inflow.  When  flux  faces  reside  on  the  surface  of  a  cutout,  the 
adjacent  out  cells  (inside  the  cutout)  may  likewise  contain  updated  values 
for  variables  entering  the  grid  from  that  location. 

At  designated  flux  boundaries,  the  velocity  component  normal  to  the 
boundary  may  be  specified  for  each  flux  face,  along  with  the  three 
Cartesian  velocity  components  for  the  out  cell  just  upstream.  In  addition, 
fixed  inflow  values  may  be  specified  for  all  the  other  transported  variables 
in  the  same  out  cell.  For  suspended  materials  (sediments  or  contaminants) 
the  concentration  in  the  out  cell  may  correspond  to  a  previous  (or  current) 
concentration  at  some  other  location  in  the  flow.  This  is  discussed  later  in 
the  section  on  reentrant  boundaries. 
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Open  boundaries  are  subject  to  an  extrapolation  procedure  developed  by 
Orlanski  (1976)  to  eliminate  spurious  reflections  in  situations  (such  as 
buoyant  flow)  that  may  generate  internal  waves.  For  normal  components 
of  velocity,  required  on  an  open  boundary,  and  for  cell-centered  variables, 
required  in  an  out  cell  just  outside  an  open  boundary,  the  Orlanski 
procedure  replaces  the  governing  equation  for  the  dependent  variable  with 


dcp  d(p  _ 

dt  dr/ 


(4.15) 


where 

c  =  effective  “wave”  speed  for  dependent  variable  (/> 

77  =  coordinate  direction  normal  to  boundary  in  computational 
space 

Equation  4.15  is  first  rearranged  to  obtain  an  expression  for  the  wave 
speed;  i.e., 


c  =  -^  (4.16) 

•A, 

Using  superscript  n  to  indicate  time  level,  and  subscript./  to  indicate 
77-position,  the  value  of  c  is  then  approximated  (for  time  level  n  -  1  and 
position/)  as 


A 77  <P--<PTX 
At  (p^-cp^ 


(4.17) 


According  to  Equation  4.15,  for  an  out  cell  at  position^'  +  1  with  an  open 
boundary  at  position^  +  1/2,  the  value  of  (f>  at  time  level  n  +  1  can  be 
approximated  as 


n+1  n  _n- 1  At  (  n  n  \  /  «  a  o\ 

<Pj+l  ~  (Pj+l~Cj  )  (4.18) 

In  the  PAR3D  version  of  the  Orlanski  procedure,  however,  this  result  is 
replaced  with 


<1  « <P"  for  c"-1  >  0 


(4.19a) 
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<P^«<PU  for  c^<  0  (4.19b) 

For  an  out  cell  at  position^  -  l  with  an  open  boundary  at  position  j  -  1/2, 
the  same  procedure  gives 


„n+i 

tPi-i 


<P> 


n+i 

Vi- 1 


‘Pi- 


tor  c"~l  <  0  (4.20a) 

for  c"^>  0  (4.20b) 


Equations  4.19  and  4.20  ensure  that  outward  propagating  disturbances 
exit  the  grid,  without  reflection,  through  the  designated  open  boundaries. 

When  used  for  cell-centered  variables,  this  extrapolation  procedure  is 
straightforward,  and  it  can  be  applied  without  further  adjustment.  When 
used  for  the  face-centered  flux  components  ( u,v,w )  defined  by 

Equations  3.22  to  3.24,  however,  an  additional  step  is  necessary. 
Specifically,  after  the  flux  extrapolation  has  been  completed,  the  resulting 
(total)  flow  rate  is  computed  for  the  contiguous  cell  faces  that  define  each 
(separately  designated)  open  boundary.  This  flow  rate  may  be  different 
from  the  specified  (fixed)  flow  rate  for  that  boundary,  and  the  ratio  of  fixed 
rate  to  extrapolated  rate  is  used  as  a  correction  factor  for  the  individual 
fluxes.  Thus,  the  distribution  of  flow  may  change  with  time,  but  the  rate  of 
flow  remains  constant  for  the  open  boundary  in  question. 

Note  that  the  flux  designation  can  be  used  for  inflow  and  outflow 
boundaries  alike,  but  the  open  designation  is  generally  limited  in  practice 
to  outflow  boundaries.  If  an  open  designation  is  to  be  used  successfully  for 
inflow,  the  boundary  in  question  should  lie  far  upstream  of  any  perturbing 
obstacles  in  the  flow. 

Reentrant  boundaries 

Inflow  boundaries  with  a  flux  designation  may  also  serve  as  reentrant 
boundaries  for  suspended  materials  (sediments  or  contaminants).  Out 
cells  adjacent  to  an  inflow  may  be  referred  to  a  single  (user-designated) 
grid  cell,  from  which  a  current  or  previous  material  concentration  is 
extracted  to  serve  as  the  updated  out-cell  concentration.  The  aim  of  this 
option  is  to  facilitate  the  modeling  of  transport  in  closed  systems  that 
recycle  suspended  materials  through  networks  of  pipes  or  ducts.  The 
reference  (intake)  grid  cell  might  lie  next  to  an  intake  for  a  pump  or  an  air 
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handler,  where  a  record  of  the  concentration  is  kept  from  one  time-step  to 
the  next.  To  account  for  travel  time  through  connecting  pipes  or  ducts,  the 
model  user  can  activate  a  delay  option,  which  computes  return  time  tR 
using  the  return  velocity  vr  and  the  return  distance  Ir  from  the  intake  cell; 
i.e., 


tR=—  (4.21) 

VR 

The  return  velocity  vr  is  assumed  to  be  the  local  inflow  velocity,  and  the 
return  distance  Ir  is  taken  to  be  the  sum  of  the  Cartesian  return-distance 
components  (xr,  yR,  zr )  from  the  intake  to  the  out  cell,  such  that 

lR=XR+yR+ZR  (4-22) 

Equation  4.22  is  intended  to  approximate  return  distances  associated  with 
connecting  pipes  and  ducts  that  follow  straight  (x,  y,  z )  lines  except  when 
changing  direction  via  right-angle  turns.  Once  the  return  time  has  been 
calculated,  the  out-cell  concentration  at  time  t  is  equated  with  the  intake 
concentration  at  time  t  -  tR. 

PAR3D  also  allows  out  cells  to  be  referred  individually  to  nearby  intake 
cells  from  which  the  suspended-material  concentration  is  transferred 
without  accounting  for  travel  time.  The  combination  of  idealized  transfers, 
with  and  without  delay,  allows  the  code  to  model  the  role  of  pipes  and 
ducts  in  closed  systems  without  computing  the  internal  details  of  the  pipe 
or  duct  flow. 

Air-water  interfaces 

For  applications  in  deep  water,  the  overhead  boundary  may  be  designated 
a  free  surface.  At  user-specified  time-intervals,  PAR3D  displaces  the  free 
surface  by  an  incremental  distance  Az  =  w  At,  where  w  is  the 
z-component  of  the  flow  velocity  normal  to  the  surface.  To  accommodate 
this  displacement,  the  z-coordinates  of  all  the  grid  nodes  are  adjusted 
proportionately  in  the  grid  block  adjacent  to  the  surface.  When  the 
computed  flow  reaches  steady  state,  the  fluid  velocity  is  tangent  to  the  free 
surface  ( w  =  o)  and  there  is  no  further  displacement  thereof. 
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The  PAR3D  treatment  of  free-surface  deformation  for  water  is  decidedly 
quasi-static.  Each  incremental  displacement  of  the  surface  represents  an 
iteration  toward  a  final  (steady  state)  configuration  in  which  the  surface  is 
perfectly  aligned  with  the  flow.  The  main  limitation  is  that  dynamic 
surface  motion  cannot  be  simulated  in  this  manner,  and  standing  waves 
are  the  only  surface  waves  that  can  be  accommodated. 

The  transported  constituents  that  may  cross  a  free  surface  are  heat, 
dissolved  gas,  and  undissolved  gas  (released  from  rising  bubbles).  The 
boundary  condition  for  heat  flux  is  the  same  as  that  given  by  Equation  4.2 
for  slip  boundaries.  The  gases  under  consideration  here  are  oxygen  and 
nitrogen,  and  these  are  transported  separately  in  distinct  (dissolved  and 
undissolved)  phases.  For  dissolved  gases  in  particular,  the  surface  flux  q 
into  the  water  is  given  by 


Q  =  KL\n\ (Xs~x)  (4-23) 

where  x  is  the  dissolved-gas  concentration,  and/s  is  the  concentration  at 
saturation.  The  transfer  coefficient  Kl  is  based  on  the  concept  of 
“penetration”  as  discussed  by  Azbel  (1981)  in  which 

(4.24) 


where 

Do  =  molecular  diffusivity  of  dissolved  gas  in  water 
co  =  exchange  frequency 

The  formula  used  for  the  free-surface  transfer  coefficient  in  PAR3D  is 
equivalent  to  one  developed  by  Fortesque  and  Pearson  (1967), 


(4.25) 


with  co  related  to  the  turbulence  energy  k  and  dissipation  rate  s  by 


co*  4.62- 
k 


(4.26) 
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Note  that  the  gas-transfer  equation  (Equation  4.23)  also  applies  at  the 
surface  of  a  bubble,  albeit  with  a  different  penetration-based  formula  for 
the  transfer  coefficient.  In  this  case,  the  exchange  frequency  is  assumed  to 
be 


co  ■ 


We 

A 


(4.27) 


where 

wb  =  bubble  rise  velocity,  relative  to  water 

The  mean  (local)  separation  distance  between  the  centers  of  the  rising 
bubbles  is 


A  =  dB 


63 


n 


\ -1/3 


(4.28) 


where 

di>>  =  effective  bubble  diameter 
3  =  local  void  ratio  for  bubbles 

Accordingly,  the  formula  for  the  bubble-surface  transfer  coefficient 
becomes 


K, 


I'Do'O 

2 

( 63 N 

V  dB  ) 

l  *  ) 

(4.29) 


Note  that  the  term  “local  void  ratio”  refers  here  to  the  ratio  of  bubble 
volume  to  fluid  volume  in  each  grid  cell. 

Solid  obstacles 

Solid  obstacles  may  be  incorporated  in  the  grid  as  internal  cutouts,  as 
shown  in  Figure  3,  or  they  may  be  included  as  empty  regions  of  Cartesian 
(x,  y,  z )  space  surrounded  by  active  grid  blocks.  Since  each  grid  block  bears 
its  own  computational  (i,j,  k)  space,  the  latter  option  avoids  the  creation 
of  unused  (cutout)  grid.  The  obstacle  surfaces  are  slip  (or  no-slip ) 
boundaries  by  default,  as  discussed  previously;  but  they  can  be 
individually  designated  as  flux  boundaries  or  as  open  boundaries. 
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When  obstacles  are  included  inside  grid  blocks  as  cutouts,  the  obstacle 
surfaces  conform  to  internal  grid  lines  of  the  blocks  in  question.  When 
included  as  empty  regions  between  blocks,  however,  their  surfaces  are 
formed  by  the  corresponding  external  surfaces  of  the  bounding  blocks.  The 
latter  option  affords  more  geometric  flexibility,  because  it  allows  the 
computational  (i,j,  k )  axes  of  each  bounding  block  to  be  oriented 
differently  with  respect  to  the  fixed  Cartesian  (x,  y,  z)  axes. 

Grid-block  interfaces 

Interfaces  between  contiguous  grid  blocks  (Figure  3)  are  treated  as  solid 
boundaries  unless  designated  otherwise  in  the  model  input.  When  a  grid  is 
divided  into  multiple  blocks  to  facilitate  parallel  processing,  however, 
block  interfaces  may  also  represent  connecting  surfaces  (“grid  cuts”) 
across  which  flow  and  transport  are  continued  without  hindrance  or 
interruption. 

Since  input  designations  for  grid  cuts  can  be  rather  lengthy  for  grids  with 
many  blocks,  PAR3D  input  (other  than  the  grid  itself)  is  divided  into  two 
input  files.  The  primary  input  file  contains  everything  except  grid-cut 
designations,  whereas  the  secondary  input  file  contains  only  grid-cut 
designations.  To  reduce  the  drudgery  associated  with  the  preparation  of 
grid-cut  input  by  hand,  a  companion  code  (SPLICE)  has  been  developed, 
which  locates  the  connecting  interfaces  between  grid  blocks,  and  creates  a 
(provisional)  secondary  input  file  in  which  all  block  interfaces  are  grid  cuts 
by  default.  This  file  can  be  edited  by  the  model  user  to  eliminate  individual 
interfaces  (or  portions  thereof)  that  are  not  intended  as  grid  cuts.  A  second 
companion  code  (GINPUT)  has  been  developed  that  creates  primary  input 
for  other  boundary  types  from  user-specified  surface  coordinates. 

Sources  and  sinks 

PAR3D  allows  user-selected  groups  of  grid  cells  to  be  designated  as 
sources  from  which  (or  into  which)  heat,  undissolved  gas,  and  suspended 
material  may  emanate  (or  vanish).  These  sources  (or  sinks)  operate  for  the 
duration  of  the  simulated  time-interval,  after  which  they  can  be 
deactivated  or  moved  to  another  location  for  a  subsequent  (hot  start)  time- 
interval.  In  the  case  of  undissolved  gas  (bubbles)  in  water,  a  free  surface 
also  acts  as  a  sink,  removing  the  gas  from  the  water  at  the  rate  that  it 
arrives  at  the  surface. 
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5  Model  Operation 

Cold  start 

Initial  (mass-conserving)  flow  is  computed  by  PAR3D  from  user-specified 
velocity  data  via  Equations  3.18  and  3.19.  Ordinarily  these  velocity  data  are 
specified  exclusively  along  inflow  and  outflow  boundaries,  in  which  case 
the  computed  initial  flow  is  potential  flow.  Any  user-specified  velocity  field 
is  admissible,  however,  and  the  resulting  initial  flow  preserves  inherent 
vorticity  while  eliminating  continuity  violations. 

Starting  with  the  initial  flow  posed  at  time  zero,  PAR3D  uses 
Equations  2.1,  2.2,  2.7,  and  2.8  (in  discrete  form)  to  advance  the 
developing  flow,  through  time,  toward  steady  state.  If  no  steady  state 
exists  for  the  situation  under  consideration,  the  computed  flow  may  reach 
a  periodic  or  quasi-periodic  condition  (e.g.,  vortex  shedding)  that  persists 
forever. 

For  passive  constituents  emanating  from  a  fixed  source,  with  null 
concentrations  initially,  transport  via  Equation  2.5  may  be  delayed  until 
the  source  itself  is  activated.  Likewise,  for  constituents  introduced  through 
inflow  boundaries,  transport  may  be  delayed  until  the  initiation  of 
introduction. 

Hot  start 

Flow  simulations  are  rarely  completed  in  a  single  run.  They  are  usually 
executed  in  a  sequence  of  runs,  called  hot  starts,  which  continue  any  flow 
development  that  was  left  incomplete  by  the  previous  run. 

Before  executing  a  hot  start,  PAR3D  reads  an  output  file  from  the  previous 
run,  which  contains  current  values  for  the  grid  coordinates  (which  may  be 
different  from  the  initial  grid)  and  all  the  dependent  variables.  It  then 
reads  the  primary  and  secondary  input  files,  which  may  direct  the  code  to 
continue  execution  with  the  same  boundary  conditions  or  with  new 
boundary  conditions.  The  latter  option  affords  users  the  freedom  to 
change  boundary  conditions  (and  source  terms)  at  will  from  one  hot  start 
to  the  next. 
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Source  code 

PAR3D  is  written  in  FORTRAN  77  with  the  MPI  Library  included  during 
compilation.  MPI  commands  embedded  in  the  source  code  mediate  the 
transfer  of  shared  information  between  processors.  Thus,  PAR3D  can  be 
executed  only  on  computers  equipped  with  FORTRAN  compilers  and  MPI 
libraries.  User-prepared  input  is  specified  in  NAMELIST  format,  and  grids 
are  accepted  in  PLOT3D  format. 
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6  Applications 

Previous  applications 

The  present  work  is  the  first  ERDC  report  documenting  PAR3D  in 
particular,  but  the  model  has  been  used  since  1999  in  a  number  of 
research  efforts.  These  include  studies  involving  secondary  flow  in  curved 
channels  (Bernard  et  al.  1999),  buoyant  flow  and  gas  transfer  induced  by 
bubble  plumes  (Bernard  2003),  and  prevention  of  salinity  intrusion  in  the 
Panama  Canal  (LaBonte  and  Spani  2005).  Since  past  work  has  been 
concerned  with  flow  and  transport  exclusively  in  water,  the  applications 
presented  here  deal  with  flow  and  transport  in  air. 

Air  flow  and  aerosol  dispersion  in  a  simple  enclosure 

Shimada  et  al.  (1996)  have  conducted  experiments  concerning  the  forced- 
air  transport  of  fine  particles  in  a  ventilated  room.  The  room  under 
consideration  (Figure  4)  was  designed  to  simulate  flow  patterns  inside  a 
train  car,  with  a  supply  vent  running  the  length  of  the  ceiling,  and  a  return 
vent  in  one  of  the  walls.  The  nominal  room  dimensions  were  270  cm  in  the 
x-direction,  330  cm  in  the  y-direction,  and  200  cm  in  the  z-direction.  The 
supply  vent  was  a  slot,  5  cm  wide,  in  the  elliptical  supply  duct;  and  the 
return  vent  was  a  square  opening,  40  cm  across,  centered  at  x  =  135  cm, 
y  =  o,  and  z  =  35  cm.  The  flow  rate  was  257  L/sec,  with  the  inflow  velocity 
profile  shown  in  Figure  4. 

Figure  5  shows  the  bounding  surfaces  for  the  grid  used  in  the  PAR3D 
representation  of  the  ventilated  room.  Figures  6  and  7  show  the  central  xz- 
plane  and  the  central  yz-plane,  respectively,  for  the  same  grid.  The  grid 
spacing  ranged  from  2.5  to  5.0  cm,  and  the  grid  was  divided  into  four 
blocks  with  cutting  planes  perpendicular  to  the  y-axis.  Figures  8  and  9 
depict  computed  (steady-state)  velocity  vectors  for  the  central  xz-plane 
and  the  central  yz-plane,  respectively.  Figures  10, 11,  and  12  show 
comparisons  of  computed  and  measured  z-components  of  velocity  in  the 
central  xz-plane. 
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Figure  4.  Outline  of  ventilated  room  used  in  experiments  conducted  by  Shimada  et  al.  (1996) 
showing  overhead  supply  vent  and  sidewall  return  vent. 


Figure  5.  Bounding  grid  surfaces  used  in  PAR3D  representation  of  ventilated  room. 


ERDC/CHL  TR-07-9 


37 


Figure  7.  Central  yz-plane  for  PAR3D  representation  of  ventilated  room. 
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Figure  8.  Computed  velocity  vectors  for  central  xz-plane  of  ventilated  room. 


Figure  9.  Computed  velocity  vectors  for  central  yz-plane  of  ventilated  room. 
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Figure  10.  Comparison  of  experimental  data  for  velocity  (•)  with  PAR3D  results  computed  at 
z  =  20  cm  in  central  xz-plane,  using  slip  condition  ( - )  and  no-slip  condition  ( — ). 


Figure  11.  Comparison  of  experimental  data  for  velocity  (•)  with  PAR3D  results  computed  at 
z  =  100  cm  in  central  xz-plane,  using  slip  condition  ( - )  and  no-slip  condition  ( — ). 
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Figure  12.  Comparison  of  experimental  data  for  velocity  (•)  with  PAR3D  results  computed  at 
z  =  140  cm  in  central  xz-plane,  using  slip  condition  ( - )  and  no-slip  condition  ( — ). 

The  PAR3D  results  shown  in  Figures  8  to  12  were  obtained  using  small 
values  for  the  turbulence  energy  (k  =  1.0  cm2/sec2)  and  eddy  viscosity 
(vt  =1.0  cm2/sec)  entering  the  room  through  the  supply  vent.  Increasing 
or  decreasing  these  values  by  a  factor  of  10  had  little  effect  on  the 
computed  flow,  which  ultimately  reached  steady  state  in  all  cases.  The 
results  shown  for  the  no-slip  condition  were  obtained  by  assuming  the 
walls,  floor,  and  ceiling  to  be  hydraulically  smooth.  There  is  little 
difference  between  the  slip  results  and  the  no-slip  results,  mainly  because 
the  influence  of  friction  is  unimportant  at  distances  far  from  the 
boundaries. 

In  the  particle-dispersion  experiments,  particle  emitters  were  placed  at  the 
intersection  of  the  central  xz-  and  yz-planes,  50  cm  above  the  floor 
(emitter  A)  and  directly  on  the  floor  (emitter  B),  as  indicated  in  Figure  13. 
In  each  experiment,  the  air  flow  in  the  room  was  allowed  to  reach  steady 
state  before  the  particle  emitter  was  activated,  and  the  emitter  was 
operated  for  600  sec  before  being  shut  off.  Particle  counters,  located 
62  cm  above  the  floor  in  the  central  xz-plane,  indicated  a  sharp  rise  in 
particle  concentrations  during  the  first  100  sec  or  so.  These  concentrations 
leveled  off  within  the  next  100  sec,  and  they  remained  (roughly)  constant 
until  the  emitter  was  shut  off. 
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Figure  13.  Emitter  locations  A  and  B  for  particle-dispersion  experiments. 


PAR3D  simulations  of  the  particle  dispersion  followed  the  same 
qualitative  pattern,  with  concentrations  increasing  only  slightly  after  the 
first  200  sec.  Figures  14  to  17  show  color  maps  of  predicted  steady-state 
concentrations  in  the  central  xz-  and  yz-planes;  and  Figures  18  and  19 
show  comparisons  of  predicted  and  measured  concentration  profiles  for 
emitters  A  and  B,  respectively.  The  concentrations  presented  here  are 
normalized  values,  as  reported  by  Shimada  et  al.  (1996),  equivalent  to  a 
delivery  rate  of  0.0167  unit  per  sec  at  the  emitter.  The  PAR3D  predictions 
were  made  using  a  null  settling  velocity,  and  a  Prandtl-Schmidt  number  of 
unity,  in  the  transport  equation  (Equation  2.5). 

PAR3D  overpredicts  observed  concentrations  for  emitter  A  by  20  to 
50  percent,  and  for  emitter  B  by  100  to  300  percent.  This  is  admittedly  a 
gross  over-prediction  of  the  data,  especially  for  emitter  B;  but  it  might  be 
remedied  somewhat  by  increasing  the  assumed  value  of  the  Prandtl- 
Schmidt  number,  and  possibly  by  including  a  settling  velocity  in  the 
transport  equation  for  the  particles.  Shimada  et  al.  assumed  their  particles 
to  be  neutrally  buoyant,  however,  and  they  made  no  measurements  of 
settling  velocity  in  still  air. 
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Figure  14.  Steady-state  particle  concentration  computed  in  central  xz-plane  for  emitter  A. 


Figure  15.  Steady-state  particle  concentration  computed  in  central  xz-plane  for  emitter  B. 
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Figure  16.  Steady-state  particle  concentration  computed  in  central  yz-plane  for  emitter  A. 


Figure  17.  Steady-state  particle  concentration  computed  in  central  yz-plane  for  emitter  B. 
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Figure  18.  Comparison  of  experimental  particle-concentration  data  (•)  for  emitter  A  with 

PAR3D  results  computed  at  z  =  62  cm  in  the  central  xz-plane,  using  slip  condition  ( - )  and 

no-slip  condition  ( — ). 


Figure  19.  Comparison  of  experimental  particle-concentration  data  (•)  for  emitter  B  with 

PAR3D  results  computed  at  z  =  62  cm  in  the  central  xz-plane,  using  slip  condition  ( - )  and 

no-slip  condition  ( — ). 


Air  flow  and  aerosol  dispersion  in  a  complex  enclosure 

Whicker  et  al.  (2000)  used  CFD  to  determine  effective  locations  for 
continuous  air  monitors  (CAMs)  in  a  ventilated  workroom  of  the  main 
plutonium  facility  at  Los  Alamos  National  Laboratory  (LANL).  The 
workroom  under  consideration  (Figure  20)  contained  five  rows  of  glove 


ERDC/CHL  TR-07-9 


45 


boxes  and  an  overhead  trolley.1  Ventilation  was  provided  by  four  inlet  air 
diffusers  in  the  ceiling  and  four  exhaust  vents  near  the  floor,  as  shown  in 
Figure  21.  The  dimensions  of  the  trolley  and  the  glove  boxes  are  shown  in 
Figure  22.  To  validate  their  CFD  code  prior  to  the  CAM  placement  study, 
Whicker  et  al.  conducted  experiments  in  which  aerosol  particles  were 
released  at  one  location  and  monitored  at  three  other  locations  in  the 
workroom.  An  overhead  view  of  the  release  and  detector  stations  is  shown 
in  Figure  23.  The  particle  source  (station  11)  was  placed  4.5  ft  above  the 
floor,  and  the  particle  detectors  (counters)  for  stations  8, 13,  and  15  were 
placed  at  heights  of  7.0,  2.0,  and  10.8  ft,  respectively. 


Figure  20.  Isometric  view  of  the  LANL  plutonium-facility  workroom  used  for  dispersion 
experiments  by  Whicker  et  al.  (2000). 


1  Figures  20  to  23,  which  show  the  experimental  workroom  layout,  are  taken  directly  from  the  LANL 
report  by  Whicker  et  al.  (2000). 
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Figure  21.  Overhead  view  of  the  workroom. 


During  the  LANL  dispersion  experiments,  the  room  ventilation  system 
was  operated  with  a  constant  (total)  flow  rate  of  approximately  85.3  cfs, 
which  was  equivalent  to  10  room-air  exchanges  per  hour.  The  inlet 
diffusers  were  circular  with  a  diameter  of  1  ft,  and  these  injected  air 
laterally  into  the  room  at  a  28-degree  angle  relative  to  the  ceiling.  The 
exhaust  vents  were  1.33  ft  wide  and  2.0  ft  high,  and  centered  2.0  ft  above 
the  floor.  The  flow  rates  through  the  four  inlet  diffusers  were  roughly  the 
same  (about  21.33  cfs  each),  and  the  average  flow  rates  through  exhaust 
vents  1,  2,  3,  and  4  were  25.6,  24.7, 14.5,  and  20.5  cfs,  respectively.1 


1  Flow  rates  were  provided  by  S.  Konecni,  co-author  of  the  report  by  Whicker  et  al.  (2000),  who  prepared 
and  executed  the  CFD  simulations  for  the  CAM  placement  study. 
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Figure  22.  Dimensions  for  glove  boxes  (GB)  in  the  workroom. 


The  computational  grid  used  for  the  PAR3D  simulations  had  an  average 
spatial  increment  of  0.5  ft  in  all  three  directions,  with  slightly  smaller 
increments  around  the  inlet  diffusers.  Figure  24  shows  overhead  views  of 
the  grid  at  different  heights  (z)  relative  to  the  floor,  and  Figure  25  shows 
side  views  of  the  grid  at  the  glove  boxes.  The  grid  was  sufficiently  fine  to 
resolve  primary  eddies  in  the  simulated  flow  field,  but  not  fine  enough  to 
resolve  boundary  layers  along  the  walls  and  the  glove  boxes.  Therefore,  all 
solid  surfaces  were  modeled  as  slip  (frictionless)  boundaries. 
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Figure  23.  Release  and  detector  stations  for  the  particle-dispersion  experiments. 


The  circular  (ceiling-mounted)  inlet  diffusers  were  represented  as 
rectangular  insets  in  the  PAR3D  grid.1  The  inflow  distributions  were 
adjusted  to  achieve  a  roughly  circular  (axially  symmetric)  flow  pattern  in 
the  radial  direction,  with  a  28-degree  angle  relative  to  the  ceiling.  The 
resulting  velocity  vectors  are  shown  around  the  base  of  Inlet  Diffuser  1  in 
Figure  26. 

Assuming  that  the  inlet  diffusers  were  supplied  by  smooth-walled,  l-ft- 
diameter,  circular  ducts,  and  that  the  supply  flows  (21.33  cfs)  were  fully 
turbulent,  then  Equations  4.7  and  4.10  imply  that  the  supply-flow 
turbulence  energy  would  have  been  k  =  3.42  ft2/sec2.  Further  assuming 
that  Equation  4.11  approximates  the  local  dissipation  rate  in  the  ducts 
(with  <5i  representing  distance  from  the  walls),  then  the  average  dissipation 


1  The  floor  is  located  atz  =  0.0  ft,  and  the  ceiling  atz  =  13.6  ft.  The  xy-locations  of  the  inlet  diffusers  are 
indicated  by  the  four  small  openings  in  the  ceiling  grid  in  Figure  24,  and  the  downward  protrusion  of 
the  modeled  diffusers  is  shown  in  Figures  25  and  26.  The  (solid)  base  of  each  diffuser  was  three  grid 
spaces  square,  with  a  total  base  area  of  1.0  ft2.  Each  diffuser  protruded  two  spaces  (0.67  ft) 
downward  into  the  grid,  and  air  was  injected  laterally  from  the  lower  cell  faces  adjacent  to  the  base. 

The  upper  lateral  faces  of  the  diffuser  were  represented  as  solid  faces,  which  allowed  return  flow  along 
the  ceiling  to  turn  downward  at  the  diffuser. 
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rate  would  have  been  e  =  15.6  ft2/sec3,  which  is  equivalent  to  an  eddy 
viscosity  vt  =  0.0675  ft2/sec. 
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Figure  24.  Overhead  views  of  computational  grid  at  different  heights  (z)  above  the  floor  of  the 

workroom. 
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Figure  25.  Side  views  of  computational  grid  (looking  west)  at  lateral  positions  (y)  relative  to 

the  east  wall  of  the  workroom. 


Figure  26.  Overhead  and  side  views  of  computed  velocity  vectors  around  the  base  of  Inlet 

Diffuser  1  for  the  workroom. 


These  values  for  supply-flow  turbulence  energy  and  eddy  viscosity 
represent  approximate  lower  bounds  for  the  workroom  inflow  values  at 
the  diffusers.  When  they  were  used  in  a  PAR3D  simulation  of  the  room  air 
flow,  however,  no  steady  state  was  achieved.  After  about  10  min  of 
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simulated  time,  the  flow  reached  a  quasi-periodic  condition  that  was  well- 
behaved  but  persistently  time-varying.  Moreover,  reducing  the  time-step 
Af  (from  its  nominal  value  of  0.02  sec)  did  not  eliminate  the  unsteadiness; 
it  merely  provided  better  time-resolution  of  the  unsteady  behavior. 

When  the  inflow  value  for  eddy  viscosity  was  increased  to  Vt  =  0.10  ft2/sec, 
the  flow  remained  unsteady;  but  when  it  was  increased  to  Vt  = 

0.15  ft2/sec,  the  flow  approached  steady  state  after  about  20  min  of 
simulated  time.  When  these  two  values  were  used  as  fixed  (uniform) 
viscosities  throughout  the  workroom,  the  lower  value  produced  unsteady 
flow,  and  the  higher  value  produced  steady  flow. 

In  retrospect,  perhaps,  it  should  not  be  surprising  that  there  existed  a 
viscosity  threshold  for  steady-state  flow  in  the  room  under  consideration. 
The  glove  boxes  and  the  overhead  trolley  all  represent  obstacles  that  could 
individually  trigger  vortex  shedding  in  a  simpler  (unidirectional)  flow 
situation;  but  the  configuration  of  the  workroom  is  so  complex  that  one 
cannot  guess  the  spatial  distribution  or  the  temporal  behavior  of  the  flow 
in  advance.  In  any  case,  increasing  the  value  of  Vt  beyond  0.15  ft2/sec 
ensured  that  PAR3D  would  yield  a  steady  state  for  the  workroom 
discretized  as  shown  in  Figures  24  and  25. 

In  the  absence  of  experimentally  determined  values  for  the  turbulence 
energy  and  the  eddy  viscosity  (or  alternatively  the  dissipation  rate)  at  the 
diffusers,  results  are  presented  that  were  obtained  using  k  =  3.42  ft2/sec2 
and  vt  =  0.15  ft2/sec  as  the  inflow  boundary  values.  Using  the  flow  rates 
outlined  above  for  the  inlet  diffusers  and  the  exhaust  vents,  PAR3D  was 
executed  for  30  min  of  simulated  time  to  ensure  that  the  computed  flow 
reached  steady  state.  Figures  27  to  30  present  overhead  views  of  the 
computed  velocity  vectors  at  four  different  heights  above  the  floor,  and 
Figure  31  offers  side  views  at  four  lateral  stations  between  the  glove  boxes. 
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Figure  27.  Overhead  view  of  computed  velocity  vectors  2.0  ft  above  floor  of  the  workroom. 


Figure  28.  Overhead  view  of  computed  velocity  vectors  4.2  ft  above  floor  of  the  workroom. 
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Figure  30.  Overhead  view  of  computed  velocity  vectors  8.3  ft  above  floor  of  the  workroom. 


ERDC/CHL  TR-07-9 


54 


y=42.1  ft 


y=28.2  ft 


y=  7.0  ft 


Figure  31.  Side  views  (looking  west)  of  computed  velocity  vectors  at  intermediate  lateral 

positions  (y)  between  glove  boxes. 


In  the  LANL  experiments,  neutrally  buoyant  aerosol  particles  were 
released  from  station  n  at  a  constant  rate  of  o.i  kg/sec  for  30  sec. 
Following  the  same  procedure  with  PAR3D,  a  point  source  was  activated 
at  station  11  after  30  min  of  flow  simulation,  which  released  contaminant 
at  0.1  kg/sec.  After  30  sec  of  simulated  release,  the  source  was  shut  off,  but 
the  dispersion  simulation  (with  an  assumed  Prandtl-Schmidt  number  of 
unity)  was  continued  for  an  additional  360  sec. 
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Figures  32,  33,  and  34  present  color  maps  of  particle  concentrations  at 
heights  of  2.0,  7.0,  and  10.8  ft,  predicted  120  sec  after  source  activation 
(and  90  sec  after  shutoff).  These  are  the  heights  of  detector  stations  13,  8, 
and  15,  respectively,  for  which  predicted  and  observed  time-histories  are 
compared  in  Figures  35,  36,  and  37.  Note  that  the  instant  represented  by 
the  color  maps  (120  sec  after  starting  release)  falls  between  the  times  at 
which  the  peak  concentrations  were  predicted  for  stations  8  and  13. 


Figure  32.  Overhead  view  of  normalized  particle  concentration,  2  ft  above  floor  of  workroom, 
120  sec  after  beginning  of  release,  predicted  by  PAR3D. 
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Figure  33.  Overhead  view  of  normalized  particle  concentration,  7  ft  above  floor  of  workroom, 
120  sec  after  beginning  of  release,  predicted  by  PAR3D. 
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Figure  34.  Overhead  view  of  normalized  particle  concentration,  10.8  ft  above  floor  of 
workroom,  120  sec  after  beginning  of  release,  predicted  by  PAR3D. 
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Figure  35.  Comparison  of  normalized  particle-concentration  data  (•)  with  PAR3D  predictions 
( — -)  for  station  13,  located  2  ft  above  workroom  floor. 


0.500 


0.400 


Figure  36.  Comparison  of  normalized  particle-concentration  data  (•)  with  PAR3D  predictions 
( — )  for  station  8,  located  7  ft  above  workroom  floor. 
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Figure  37.  Comparison  of  normalized  particle-concentration  data  (•)  with  PAR3D  predictions 
( — )  for  station  15,  located  10.8  ft  above  workroom  floor. 


PAR3D  computes  contaminant  concentration  in  units  of  mass  per  unit 
volume,  but  the  experimental  concentrations  in  the  report  by  Whicker 
et  al.  (2000)  were  given  as  particles  per  unit  volume.  No  masses  were 
given  for  the  aerosol  particles,  so  it  is  not  possible  here  to  compare  the 
computed  and  measured  concentrations  in  the  same  dimensional  units. 
The  normalized  PAR3D  predictions  in  Figures  32  to  37  were  obtained  by 
dividing  all  the  (dimensional)  predicted  concentrations  by  the  maximum 
(dimensional)  concentration  predicted  for  detector  station  13.  Likewise, 
the  normalized  experimental  data  were  obtained  by  dividing  all  the 
(dimensional)  measured  concentrations  by  the  maximum  (dimensional) 
concentration  measured  at  station  13. 

The  most  favorable  comparison  with  the  experimental  data  was  achieved 
for  station  13,  as  indicated  by  Figure  35.  The  arrival  time  for  the  peak 
concentration  is  late  by  about  20  sec,  but  the  shape  of  the  curve  compares 
well  with  the  data.  Note  that  the  predicted  peak  matches  the  observed 
peak  for  this  station  because  of  the  normalization. 

Similar  comments  apply  for  station  8,  as  demonstrated  by  Figure  36.  The 
initial  rise  in  concentration  is  late  by  about  40  sec,  but  the  peak 
concentration  is  late  by  no  more  than  20  sec.  The  shape  of  the  predicted 
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curve  also  compares  well  with  the  data,  and  the  predicted  peak  matches 
the  observed  peak  within  to  percent.  Note  that  the  agreement  between 
predicted  and  observed  peaks  for  this  station  indicates  that  the  relative 
peaks  are  correctly  predicted  for  detectors  8  and  13,  with  or  without 
normalization. 

Predicted  results  for  station  15  are  poor  in  comparison  with  those  for 
stations  8  and  13,  as  is  obvious  from  Figure  37.  The  initial  rise  in 
concentration  is  late  by  at  least  too  sec,  and  the  highest  predicted 
concentration,  which  occurs  at  the  end  of  the  simulated  time-interval,  is 
lower  than  the  observed  peak  by  a  factor  of  two  or  more. 

The  predictions  in  Figures  35  to  37  represent  the  best  that  PAR3D  was 
able  to  achieve  (with  steady  flow)  for  the  workroom  under  consideration. 

In  subsequent  simulations,  using  k  =  3.42  ft2/sec2  with  vt  >  0.15  ft2/sec  at 
the  diffusers,  it  was  found  that  the  agreement  between  prediction  and 
experiment  worsened  with  increasing  vt.  In  this  case,  the  provisional  rule 
of  thumb  (for  predicting  dispersion)  seems  to  be  that  one  should  use  an 
inflow  value  for  vt  roughly  equal  to  the  smallest  value  that  produces  steady 
flow. 
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7  Conclusions 

The  PAR3D  code,  initially  developed  as  a  numerical  model  for  simulating 
flow  and  transport  in  deep  water,  has  been  extended  for  applications 
involving  air  flowthrough  complex  buildings  and  other  enclosures. 
Comparisons  with  experimental  results  here  and  elsewhere 
(Bernard  et  al.1999;  Bernard  2003)  indicate  that  velocity  predictions  may 
be  acceptable,  but  that  dispersion  predictions  may  be  somewhat  imprecise. 
Nevertheless,  experience  thus  far  suggests  that  the  model  is  qualitatively 
correct  in  its  predictions  (most  importantly  those  presented  herein  for  a 
complex  enclosure)  and  that  it  may  be  useful  for  gaining  insight  into  the 
dispersion  of  contaminants  under  diverse  flow  conditions. 

When  the  code  is  executed  with  the  turbulence  model,  inflow  values  for 
the  turbulence  energy  and  its  dissipation  rate  (or  associated  eddy 
viscosity)  may  serve  as  free  parameters  for  quantitatively  tuning  predicted 
results.  These  quantities  can  be  determined  in  principle  from  experiments 
or  from  empirical  considerations,  but  their  values  in  practice  may  be 
uncertain  by  an  order  of  magnitude  or  more.  In  any  case,  the  main  reasons 
for  including  the  turbulence  model  are  to  capture  the  qualitative  effects  of 
turbulence  on  the  mean  flow  and  to  avoid  having  the  model  user  specify  a 
viscosity  for  the  entire  flow  field. 

When  the  code  is  executed  without  the  turbulence  model,  then  the  model 
user  must  specify  a  single  (uniform)  value  for  viscosity  that  will  be  used  for 
the  duration  of  the  simulated  time-interval.  Unless  the  aim  is  to  reproduce 
truly  laminar  flow  at  a  given  Reynolds  number,  this  viscosity  serves  as  a 
free  parameter  that  can  be  used  to  adjust  predicted  results  for  agreement 
with  experimental  data  and/or  sponsor  expectations. 

Although  free  parameters  can  be  used  (within  reasonable  bounds)  to  tune 
quantitative  predictions  made  by  PAR3D,  qualitative  verisimilitude  is 
mainly  a  matter  of  adequate  grid  resolution  in  three  dimensions.  At  one 
extreme,  if  the  grid  is  too  coarse  to  capture  the  primary  eddies  in  a  given 
setting,  then  the  resulting  flow  simulation  is  qualitatively  worthless.  At  the 
other  extreme,  grid  independence  is  achieved  when  further  refinement  of 
the  grid  produces  no  further  change  in  quantitative  predictions.  The 
degree  of  grid  refinement  employed  for  the  applications  in  Chapter  6  was 
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more  than  adequate  for  the  primary  flow  (and  related  transport)  but  not 
quite  sufficient  for  true  grid  independence. 

The  time  required  for  PAR3D  execution  depends  on  the  number  of 
processors  used,  the  (floating-point)  computational  speed  of  the  individual 
processors,  and  the  time  required  for  communication  between  processors. 
Computational  speed  and  communication  time  vary  considerably  from 
computer  to  computer,  but  they  continue  to  improve  dramatically  with 
each  new  generation  of  computers.  For  the  applications  presented  in 
Chapter  6,  the  simple  enclosure  required  several  hours  (with  four 
processors)  and  the  complex  enclosure  required  an  overnight  run  (with 
eight  processors)  on  an  SGI  Origin  3900  computer. 
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Appendix:  Boundary  Conditions  for 
Turbulence 

The  k-e  turbulence  model  (Chapter  2)  entails  only  the  governing 
equations  for  the  turbulence  energy  k  and  the  turbulence-energy 
dissipation  rate  e.  Thus,  the  boundary  conditions  for  these  two  quantities 
have  to  be  determined  from  other  considerations. 

Inflow  boundaries  are  straightforward  in  principle  because  the  inflow 
values  for  k  and  e  should  be  the  same  as  those  in  the  flow  upstream  of  the 
inflow.  In  practice,  however,  the  flow  upstream  may  be  poorly 
characterized;  this  can  render  the  inflow  conditions  somewhat 
problematic,  as  exemplified  by  the  model  applications  discussed  in 
Chapter  6.  Outflow  boundaries  present  less  of  a  problem  because  they 
usually  lend  themselves  to  the  extrapolation  procedure  (Equations  4.15  to 
4.20)  outlined  in  Chapter  4.  Slip  boundaries  are  even  less  troublesome 
because  they  can  be  idealized  as  symmetry  boundaries  through  which 
there  is  no  gradient  flux  of  k  or  e,  as  stipulated  by  Equation  4.1. 

No-slip  boundaries  (friction  walls)  present  a  different  problem  altogether 
because  they  require  in  principle  a  detailed  characterization  of  the 
turbulent  flow  near  the  walls.  This  often  demands  much  finer  grid 
resolution  than  is  feasible  in  practice,  and  one  is  forced  to  seek  an  easier 
way  out.  The  most  common  remedy  seems  to  be  the  use  of  wall  functions 
for  the  values  of  k  and  e  in  the  flow  near  the  walls,  as  discussed  by  Rodi 
(1980).  For  unidirectional  shear  flow  without  a  pressure  gradient,  these 
are  given  by 


k  = 


for  z  <  S± 


(A.l) 


s  =  —  for  z<81 

KZ 


(A. 2) 


where 

v *  =  friction  velocity 
Cv  =  0.09 
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z  =  normal  distance  from  wall 

<5i  =  thickness  of  wall-adjacent  grid  cell  normal  to  wall 
k  =  0.40 

Since  these  wall  functions  are  imposed  in  the  flow  rather  than  on  the  walls, 
they  do  not  represent  true  boundary  conditions.  In  fact,  they  are  truly 
valid  only  for  turbulent  shear  flow  in  equilibrium,  such  that 

P  =  s  (A. 3) 

where  P  is  the  turbulence-energy  production  rate,  given  by  Equation  2.11 
in  Chapter  2. 

As  an  alternative  to  wall  functions,  one  seeks  boundary  conditions  that 
yield  the  near-wall  values  for  k  and  e  (Equations  A.i  and  A.  2)  only  at 
equilibrium,  and  not  otherwise.  For  a  shear  flow  in  equilibrium,  the  eddy 
viscosity  very  near  the  wall  is  given  by 

vT=Koi,z  (A. 4) 


and  the  gradient  flux  (per  unit  area)  of  k  normal  to  the  wall  reduces  to 


vT  8k  _  k  v*  z  8k 
<jk  8 z  <jk  8z 


(A. 5) 


Assuming  that 


8  k  k  -  k* 

dz  z 


(A. 6) 


with 


k*  = 


(A. 7) 


then  Equation  A.  5  gives  way  to 


vT  8k 
crk  dz 


KV* 


(k-k*) 


(A. 8) 
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When  imposed  as  the  wall-normal  gradient  flux  at  the  wall,  Equation  A.8 
ensures  that  — >  fc*  as  the  flow  approaches  equilibrium  and  the  governing 
equation  (Equation  2.7)  reduces  to  the  limiting  condition, 


d  (  vT  dk ' 
dz{crk  dz j 


— >  0 


(A. 9) 


The  value  of  k  is  uniform  when  the  flow  is  in  equilibrium,  as  dictated  by 
Equation  A.i,  and  the  appropriate  boundary  value  is  unambiguous.  In 
contrast,  however,  Equation  A.  2  indicates  that  e  increases  sharply  as 
z  -»  o,  and  it  is  not  clear  what  boundary  value  should  be  used  to  elicit  the 
z-dependence  evident  in  the  wall  function  (A. 2).  Therefore,  one  must 
begin  by  examining  the  equilibrium  form  of  the  governing  equation 
(Equation  2.8)  for  e: 


s1  (  x  1  d 
k  cr  dz 


f  <9  A 

v 


V 


dz 


=  0 


(A. 10) 


First,  consider  the  average  dissipation  rate  in  a  wall-adjacent  grid  cell, 
approximated  with  the  expression, 


s  ~ 


P  of 

k81 


(A. 11) 


in  which  ft  >  1.  Substituting  Equation  A.i  for  k  and  Equation  A.11  for  e,  the 
first  term  in  Equation  A.  10  can  be  replaced  with 


e2  ,  x 

c 2): 


P2  C 


1/2  ot 


K 


K 


(ci  ^2 ) 


(A. 12) 


Next,  using  central  differencing  to  evaluate  the  second  term  in 
Equation  A.10,  one  obtains 


1  8 

(  ds^ 

1/  _ 

1 

'(  a  A 

1/  _ 

z=S1 

f  aA 

1/  _ 

z-0 

ae  dz 

l  7  dzj 

A 

VT  ^ 

l  8zt 

VT  3 

l  8zj 

(A. 13) 


Differentiating  Equation  A.  2  with  respect  to  z  gives 
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ds  _  of 
dz  kz2 


(A. 14) 


Multiplying  Equation  A.  14  by  Equation  A.4,  and  evaluating  the  result  at 
z  =  81  yields 


f 

VT 

V 


ds) 


dz 


y  z=s1 


with  which  Equation  A.  13  reduces  to 


(A. 15) 


1  a  1 

(  ds) 

1/  _ 

,  __  1 

of 

_ |_ 

r  ds 

1/  _ 

<*e  dz 

l  Tdz) 

'  °eSl 

*1 

VT  ~ 

l  5z 

(A. 16) 


After  substituting  Equations  A.  12  and  A.  16  for  the  first  and  second  terms, 
respectively,  in  Equation  A.10,  and  then  rearranging,  one  finds  that 


ds 

dz 


p2c 


1/2u 4 


Jz= 0 


K 


sf 


(ci  ^2 ) 


of 

Oe5l 


(A.  17) 


Next,  assuming  that 


ds  s-s* 
dz  z 


(A. 18) 


and  substituting  Equation  A.4  for  vt,  the  left  side  of  Equation  A.  17  takes 
the  form 


1 

f  ds ) 

1/  _ 

~  KO* 

f  pvt  ) 

c* 

)T  dzj 

2=0 

l  *A  j 

Combining  Equations  A.17  and  A.19,  it  follows  that 

P  of  of 

£■*  « - 1 - 

kS1  k81 


^A(c2-ci)+1 

K 


(A.19) 


(A. 20) 


Lastly,  the  equilibrium  relation,  Equation  2.9,  dictates  that 
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C2  C± 


K 


a  C 

£  V 


1/2 


(A. 21) 


so  that  Equation  A. 20  gives  way  to 

s^-^-[p2+P  +  l\  (A. 22) 

K  O ^ 

Equation  A.22  defines  the  boundary  value  e*  that  should  be  used  for 
computing  the  gradient  flux  of  e  normal  to  the  wall  in  Equations  A.  19 
and  4.4. 

What  now  remains  to  be  determined  is  the  influence  of  the  wall  on  the 
computed  rate  of  turbulence-energy  production  P.  Invoking  Equations  A.3 
and  A.11,  the  average  value  in  a  wall-adjacent  grid  cell  is  given 
approximately  by 


P» 


plK 

kS1 


(A. 23) 


and  Equations  2.6  and  A.11  imply  that  the  average  eddy  viscosity  is 


ko„81 

~T~ 


(A. 24) 


For  the  shear  flow  under  consideration,  however,  Equation  2.11  reduces  to 


P0  =  VT 


v5zy 


(A. 25) 


Due  to  the  proximity  of  the  wall,  and  the  placement  of  computed  velocities 
in  the  flow,  PAR3D  would  be  forced  to  evaluate  the  z-derivative  in 
Equation  A.25  at  the  top  of  the  wall-adjacent  grid  cell  (i.e.,  at  z  =  §1). 
Assuming  a  logarithmic  velocity  profile  (similar  to  that  in  Equation  4.7c) 
with  the  form 


u  =  — In  z  +  const 

K 


(A. 26) 
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then  the  exact  value  of  the  z-derivative  at  z  =  <5i  is 


f  du^ 

\dzJ,-Sl 


v, 


Combining  Equations  A.24,  A.25,  and  A.27  now  yields 


(A. 27) 


Po  = 


Pk8± 


(A. 28) 


which  is  the  production  rate  that  PAR3D  would  predict  for  the  wall- 
adjacent  grid  cell  (using  computed  velocities  in  the  flow  field)  if  it  could  do 
so  with  infinite  precision.  Note,  however,  that  the  production  rate  given  by 
Equation  A.  2  8  does  not  match  the  cell-averaged  production  rate  given  by 
Equation  A.23.  To  make  up  the  deficit,  Equation  A.28  (and  likewise 
Equation  2.11)  has  to  be  supplemented  with  an  additional  wall -induced 
production  term, 


Pi  = 


o* 


k81 


J3-- 

P. 


(A. 29) 


and  this  is  the  basis  for  the  additional  term  in  the  wall-adjusted 
production  rate  given  by  Equation  4.12. 


The  free  parameter  f  is  dependent  on  the  numerical  discretization  of  the 
governing  equations,  and  its  value  has  to  be  determined  (by  trial  and 
error)  from  computational  experiments  with  unidirectional,  wall-bound, 
turbulent  shear  flow,  for  which  Equations  A.i  and  A.  2  define  the 
equilibrium  values  of  k  and  e,  respectively.  For  this  type  of  flow,  in  the 
absence  of  a  pressure  gradient,  the  Reynolds-averaged  momentum 
equation  (Equation  2.1)  reduces  to 


d_ 

dz 


r  du ^ 
'T~dz 


=  0 


(A. 30) 


and  the  steady-state  governing  equations  for  the  turbulence  quantities 
become 
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„  1  5 

P~£~{ - 

ak  dz 


dk 

dz 


=  0 


(A. 31) 


s 

k 


£) 


i  a  ( 


+  a  E  dz 


d£ A 

dZy 


=  0 


(A. 32) 


Unpublished  work  by  one  of  the  authors  (Bernard),  in  which 
Equations  A.30  to  A.32  were  solved  numerically  with  a  l-D  MAC 
discretization  (Chapter  3),  indicates  that  the  correct  profiles  for  velocity, 
energy,  and  dissipation  rate  are  consistently  reproduced  within  a  few 
percent  when  the  value  used  for  the  free  parameter  is  (3  *  2.  Inserting  this 
value  for  / 3  in  Equations  A.22  and  A.29,  respectively,  gives 


(A. 33) 


O  3 

P1=-  — 
2  kS1 


(A. 34) 


The  developments  presented  above  are  intended  to  remove  any  mystery 
associated  with  the  wall-adjusted  production  rate,  given  by  Equation  4.12, 
and  the  wall-normal  gradient  fluxes  for  k  and  e,  given  by  Equations  4.13 
and  4.14,  respectively,  in  Chapter  4.  Note  that  the  questionable 
assumptions  in  this  process,  even  for  turbulent  shear  flow  in  equilibrium, 
are  those  embodied  by  Equations  A. 6  and  A.  18.  The  end  results,  however, 
are  boundary  conditions  that  adequately  reproduce  wall-bound  shear  flow 
(when  used  with  the  k-e  turbulence  model)  while  accommodating  other 
types  of  flow,  without  imposing  wall-function  values  in  the  flow  field  away 
from  the  walls.  The  latter  attribute  is  important  when  the  grid  is  too  coarse 
to  resolve  near-wall  flow  because  it  prevents  the  computed  flow  from  being 
(erroneously)  over-influenced  by  the  walls. 
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